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control  effort  results.  From  the  control  effort,  useful  design  parameters  such  as  AV  and 
propellant  mass  are  determined.  For  formation  maneuvering,  the  drone  spacecraft  track 
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ABSTRACT 


A  growing  interest  in  formation  flying  satellites  demands  development  and 
analysis  of  control  and  estimation  algorithms  for  station-keeping  and  formation 
maneuvering.  This  thesis  discusses  the  development  of  a  discrete  linear-quadratic- 
regulator  control  algorithm  for  formations  in  the  vicinity  of  the  L2  sun-earth  libration 
point.  The  development  of  an  appropriate  Kalman  filter  is  included  as  well.  Simulations 
are  created  for  the  analysis  of  the  station-keeping  and  various  formation  maneuvers  of  the 
Stellar  Imager  mission.  The  simulations  provide  tracking  error,  estimation  error,  and 
control  effort  results.  From  the  control  effort,  useful  design  parameters  such  as  AV  and 
propellant  mass  are  determined.  For  formation  maneuvering,  the  drone  spacecraft  track 
to  within  4  meters  of  their  desired  position  and  within  1.5  millimeters  per  second  of  their 
desired  zero  velocity.  The  filter,  with  few  exceptions,  keeps  the  estimation  errors  within 
their  three-sigma  values.  Without  noise,  the  controller  performs  extremely  well,  with  the 
drones  tracking  to  within  several  micrometers.  Each  drone  uses  around  1  to  2  grams  of 
propellant  per  maneuver,  depending  on  the  circumstances. 
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1.  INTRODUCTION 


1.1  Libration  Point  Introduction 

Discovered  by  Euler  and  Lagrange  in  1772,  while  studying  the  motion  of  the 
moon,  libration  points  are  a  significant  niche  in  the  field  of  orbital  dynamics.  In  the 
simplest  sense,  libration  points  are  points  in  space  where  the  gravitational  and  centrifugal 
forces  cancel.  Szebehely1  studied  the  mathematics  of  the  generic  restricted  three-body 
problem,  and  formulated  the  positions  and  behaviors  of  the  libration  points  and  the 
dynamics  of  objects  in  their  vicinity.  Barden2,  and  later  Howell3,  expanded  the  study  of 
libration  point  dynamics  to  their  transitions  and  trajectory  design. 

Farquhar4  initially  examined  control  techniques  for  the  station-keeping  of 
satellites  orbiting  a  libration  point.  Gomez,  Masdemont,  and  Simo5  have  extensively 
studied  dynamics  and  control  problems  involving  libration  points.  Scheeres  and  Vinh6 
have  researched  the  dynamics  and  control  of  relative  motion  of  two  spacecraft  in  unstable 
orbits.  Hoffman7  used  modem  control  techniques  to  perform  station  keeping  around  the 
earth-moon  collinear  libration  points. 

Many  future  space  missions  are  planning  to  use  libration  points.  In  fact,  the 
Microwave  Anisotropy  Probe8  (MAP)  began  orbiting  the  sun-earth  L2  point  on  1  October 
2001.  Other  future  libration  point  missions  include  Stellar  Imager9,  Micro  Arcsecond  X- 
Ray  Interferometry  Mission10  (MAXIM),  MAXIM  Pathfinder11,  and  possibly  Terrestrial 
Planet  Finder12.  Further  understanding  of  the  dynamics  and  control  of  satellites  in  the 
vicinity  of  libration  points  is  essential  for  the  success  of  these  and  many  other  future 
space  missions. 


l 


1.2  Formation  Flying  Introduction 

Currently,  formation  flying  spacecraft  control  is  being  extensively  researched. 
Flying  satellites  relative  to  one  another  offers  mission  possibilities  otherwise  infeasible. 
Multiple  satellites  working  together  allows  for  more  powerful  telescopes,  interferometers, 
and  other  scientific  data  collecting  instruments.  Communications  technology  may 
improve  as  well.  With  multiple  satellites,  risk  is  reduced  through  redundancy,  so  a 
failure  of  one  element  does  not  compromise  the  entire  mission.  Also,  if  the  spacecraft  are 
similar  or  identical,  assembly  line  production  can  reduce  costs. 

Formation  flying  control  can  be  performed  in  two  ways— centralized  or 
decentralized.  With  centralized  control,  one  spacecraft  or  processor  calculates  and 
commands  the  motion  of  the  entire  formation.  With  decentralized  control,  each 
spacecraft,  with  input  from  the  rest  of  the  formation,  processes  its  own  control 
requirements.  Although  this  thesis  presents  a  centralized  control  approach,  the 
decentralized  method  is  the  basis  for  many  references.  Speyer13  first  introduced  a 
decentralized  linear-quadratic-Gaussian  control  method.  Carpenter14, 15  applied  this  work 
to  formation  flying  satellites,  and  further  expanded  it  to  deal  with  both  time-invariant  and 
time-varying  systems.  Speyer’s13  method  produces  identical  results  to  the  centralized 
linear-quadratic-Gaussian  control  method,  and  it  also  minimizes  data  transmission. 

Formation  flying  satellites  in  Earth  orbit  have  been  demonstrated.  Carpenter, 
Folta,  and  Quinn16  derived  a  decentralized  framework  for  the  applicability  of  autonomous 
formation  flying  control  for  the  EO-1  mission  to  follow  Landsat-7  in  low-Earth  orbit. 
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The  next  step  is  extending  formation  flying  missions  to  other  points  in  space,  such  as  the 
sun-earth  libration  points. 

NASA  has  several  distributed  spacecraft  missions  planned  for  the  next  decade  and 
beyond.  The  Magnetospheric  Constellation17  will  study  the  magnetotail  of  the  earth. 
Stellar  Imager9  and  MAXIM10  will  image  stars  and  black  holes  respectively.  Also,  the 
Laser  Interferometer  Space  Antenna  will  detect  gravitational  waves  in  an  attempt  to 
prove  elements  of  general  relativity.  All  of  these  missions,  and  many  more,  rely  on 
formation  flying. 

1.3  Overview 

In  Section  2,  the  equations  of  motion  for  the  circular  restricted-three  body 
problem  are  developed  and  explained.  From  these  equations  of  motion,  the  libration 
points  are  discovered  and  addressed.  The  problem  is  simplified  by  linearizing  the 
equations  of  motion  about  the  collinear  libration  points.  Then,  the  orbits  unique  to 
libration  points  are  examined. 

In  Section  3,  the  open-loop  dynamics  of  a  system  are  briefly  explained.  State 
feedback  control  is  applied  to  “close  the  loop.”  The  optimal  control  is  found  using  the 
linear-quadratic-regulator  method.  Then,  the  problem  is  simplified  with  an  infinite 
horizon  assumption.  The  stability  of  the  linear-quadratic-regulator  is  proven.  In  Section 
4,  the  continuous  system  is  sampled,  and  the  discrete  optimal  control  law  is  found.  These 
concepts  are  then  applied  to  satellites  flying  in  formation  in  Section  5. 

In  Section  6,  an  observer  is  introduced  to  estimate  the  states  of  a  system  by  using 
incoming  measurements.  Augmenting  the  system  state-space  for  multiple  satellites  is 
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shown.  The  optimal  method  of  updating  the  estimates  is  given,  and  the  Kalman  filter  is 

derived.  A  common  Kalman  filter  algorithm  is  shown.  The  Kalman  filter  is  adapted  to 

accommodate  nonlinear  measurements,  with  range,  azimuth,  and  elevation  given  as  an 
example. 

The  background  to  the  Stellar  Imager  mission  is  explored  in  Section  7,  yielding 
performance  and  design  requirements.  Special  emphasis  is  given  to  the  preliminary 
formation,  orbit,  and  control  designs.  In  Section  8,  the  reference  orbit  and  dynamics 
specific  to  the  sun-earth  L2  point  are  determined  from  the  circular  restricted  three-body; 
however,  a  more  realistic  orbit  and  associated  dynamics  are  simulated  instead. 
Simulations  are  developed  to  maintain  a  satellite  on  this  desired  reference  orbit,  slew  the 
formation  to  aim  at  another  star,  and  reconfigure  the  formation  to  achieve  better  imaging 

results.  The  simulations  produce  useful  tracking  error,  estimation  error,  and  control 
effort  results. 
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2.  CIRCULAR  RESTRICTED  THREE-BODY  PROBLEM 


2.1  Equations  of  Motion 

The  classic  restricted  three-body  problem  refers  to  the  motion  of  an 
infinitesimally  small  body  in  the  presence  of  two  larger  bodies  orbiting  their  barycenter, 
or  common  center  of  mass.  Two  examples  of  the  restricted  three-body  problem  are  the 
earth-moon  and  sun-earth  systems.  For  the  earth-moon  system,  the  earth  is  the  first  body, 
the  moon,  orbiting  the  earth,  is  the  second  body,  and  a  satellite  would  be  the  third.  For 
the  sun-earth  system,  the  sun  is  the  first  body,  the  earth,  orbiting  the  sun,  is  the  second 
body,  and  a  satellite  would  be  third.  For  preliminary  analysis,  the  orbits  of  the  moon 
around  the  earth  and  the  earth  around  the  sun  are  assumed  to  be  circular,  have  a  constant 
angular  velocity,  and  have  no  variation  in  inclination  from  the  orbital  plane.  When  these 
assumptions  are  made,  the  restricted  three-body  problem  becomes  the  more  specific 
circular  restricted  three-body  problem.  This  thesis  will  focus  on  the  sun-earth  circular 
restricted  three-body  problem,  although  the  development  of  equations  will  hold  for  any 
three  bodies. 

For  the  development  of  the  equations  of  motion,  the  rotating  coordinate  system 
shown  in  Figure  1  will  be  used.  The  origin  is  at  the  barycenter  of  the  two  large  bodies. 
The  first  basis  vector,  a, ,  runs  along  the  line  between  the  two  large  bodies  in  the 
direction  of  the  smaller  body.  The  third  basis  vector,  a3 ,  is  in  the  direction  of  the  orbit 
normal.  The  cross  product  of  a3  and  a,  yields  the  second  basis  vector,  a2 .  These 
vectors  form  an  orthogonal  coordinate  system. 
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m3(X,y,Z) 


a2  | 


The  position  vector  of  the  spacecraft  in  the  rotating  frame  is 

r  =  Xal+Ya2+Za3.  (2.1) 

Writing  the  vectors  from  the  large  bodies  to  the  spacecraft  in  the  rotating  frame  gives 

/i  =  (X  +  D,  )a1  +  Ya2  +  Za3  (2.2) 

and 

r2  =(X  -  D2  )ax  +  Ya2  +  Za3 ,  (2.3) 

where  Dl  is  the  distance  from  the  barycenter  to  the  largest  body  and  D2  is  the  distance 
from  the  barycenter  to  the  smaller  body.  Note  that 

D  =  D1+D2.  (2.4) 

An  inertial  coordinate  system  is  defined  such  that  the  origin  is  the  barycenter,  but 
the  basis  vectors  do  not  rotate  with  the  system.  The  inertial  basis  vectors  n,  and  h2  lie  in 
the  same  plane  as  the  rotating  basis  vectors  a,  and  a2 ;  thus,  n3  is  identical  to  a3 .  Once 
every  revolution  n,  is  aligned  with  a, ,  and  n2  is  aligned  with  a2 .  With  respect  to  the 
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inertial  coordinate  frame,  the  two  large  bodies  are  rotating  about  their  barycenter  with  a 

constant  angular  velocity 

<u  =  [O  0  (of , 

(2.5) 

where 

(2.6) 

Therefore,  deriving  with  respect  to  the  inertial  frame 

dj  -(Oa2 

(2.7) 

a2  =-coai 

(2.8) 

d3=  0. 

(2.9) 

Taking  the  first  and  second  time  derivatives  of  the  position  vector  (Equation  2.1)  with 
respect  to  the  inertial  frame, 

r  =  (X-  (oY)ax  +  (coX  +  Y)a2  +  Za3  (2. 10) 

and 

r  =  (X  —  2jcoY  ~(o2X)a,  +  (Y  +  2wX  -eo2Y)a2  +  Za3 .  (2.11) 

Newton’s  Law  of  Universal  Gravitation  for  n  bodies  states 


xh  mi 

Fg  =  — -rj' , 


/=!  rH 


(2.12) 


where  F  is  the  total  gravitational  force  acting  on  the  ith  body,  G  is  the  universal 


gravitation  constant,  m  represents  the  masses  of  the  respective  bodies,  and  rJt  is  the 

position  vector  from  the  fh  body  to  the  ith  body.  Newton’s  Second  Law  defines  force  as 
the  time  derivative  of  momentum 
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„  dv,  dm, 
F  =  mi — — + V; — 

'  dt  '  dt 


(2.13) 


where 


(2.14) 


dvt  _  d  rt 
~d^~~dtr 


(2.15) 


The  derivatives  in  Equations  2.14  and  2.15  are  with  respect  to  an  inertial  frame. 
Assuming  constant  mass,  then 


—  =  0. 


(2.16) 


With  gravity  as  the  only  force  (the  restricted  part  of  the  restricted  three-body  problem), 
substituting  Equations  2.13  and  2.16  into  Equation  2.12,  and  dividing  by  m, ,  yields 


mi 

r>  ]  i TrJ‘ 


(2.17) 


J=\  ra 

j*i  Jl 


Expanding  for  the  three-body  problem  and  dropping  the  subscript  i  gives  the  equation  of 
motion  for  a  satellite  in  the  restricted  three-body  problem 


Gmj  Gm2 
r  =  -Tirri-T_irr2' 


(2.18) 


Substituting  Equations  2.11, 2.2,  and  2.3  into  Equation  2.18,  the  equations  of  motion  can 
now  be  written  as 


X-2(oY-(d2X=-^-^-^  Dl) 


(2.19) 
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(2.20) 


Y  +  2coX-cq2Y  =  - 


h2y 


\r2 


z 


(2.21) 


where 


flj  =  Gntj .  (2.22) 

Szebehely1  shows  that  the  centrifugal  plus  gravitational  force  potential,  U,  exists  such 
that 


and 


U  =-0)2(X2  +Y2)  +  ^-  + 

2  |n 


(2.23) 


X-2 cdY  = 


dU 

dX 


(2.24) 


Y  +  2coX  = 


dU 

dY 


(2.25) 

(2.26) 


2.2  Libration  Points 

By  setting  all  time  derivatives  in  the  equations  of  motion  to  zero,  five  libration,  or 
Lagrange,  points  can  be  calculated.  The  location  of  these  points  depend  on  the  masses 
and  distances  of  the  bodies;  however,  three  points  are  always  collinear  with  the  two  large 
bodies  (LI,  L2,  and  L3),  and  the  other  two  points  form  equilateral  triangles  with  the  two 
large  bodies  (L4,  and  L5).  Szebehely1  calls  the  point  between  the  two  masses  L2  and  the 
point  opposite  the  smaller  body  LI.  However,  most  of  today’s  literature,  including  Wie19 


9 


and  Farquhar4,  have  just  the  opposite,  with  LI  being  between  the  two  masses  and  L2 
opposite  the  smaller  body,  as  shown  in  Figure  2. 


For  the  sun-earth  system,  it  is  common  to  use  the  earth-moon  barycenter  as  the 
second  body  rather  than  the  earth  itself,  and  is  done  so  by  Szebehely1.  Hoffman7  shows 
that  the  motion  near  the  equilateral  libration  points  is  stable  when  the  ratio  of  the  smaller 
body  mass  to  the  total  system  mass  is  less  than  0.0385,  and  that  the  motion  near  the 
collinear  libration  points  is  always  unstable.  Barring  any  other  forces,  an  object’s  orbit 
around  a  stable  equilateral  libration  point  will  eventually  decay  until  it  rests  at  the 
libration  point  itself.  Szebehely1  mentions  the  Trojan  asteroids  sit  at  a  sun-Jupiter 
equilateral  libration  point.  An  object  in  an  orbit  around  a  collinear  libration  point,  on  the 
other  hand,  will  eventually  fly  away  from  the  vicinity  of  the  libration  point. 

2.3  Linearized  Equations  of  Motion  Near  Collinear  Libration  Points 

Motion  around  a  collinear  libration  point  is  more  easily  expressed  in  a  local 
coordinate  frame,  shown  in  Figure  3,  with  the  origin  located  at  the  libration  point.  The 
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radial  direction,  x,  is  collinear  with  a, ,  going  from  the  sun  (or  largest  body)  through  the 
libration  point.  The  cross-track  direction,  z,  is  parallel  and  in  the  same  direction  as  , 

the  orbit  normal.  Crossing  z  with  x  gives  the  in-track  (or  along-track)  direction,  y.  When 
circular  motion  is  assumed,  the  y  direction  is  the  tangential  velocity  direction  of  the 
libration  point  around  the  system  barycenter.  These  axes  form  an  orthogonal  coordinate 
system.  Motion  in  the  x-y  plane  will  be  referred  to  as  in-plane,  and  any  motion  in  the  z 
direction  will  be  referred  to  as  out-of-plane. 


Figure  3.  Relative  position  with  respect  to  L2 
Next,  break  the  position  vector,  r,  into  the  libration  point  position  vector,  r0,  and 
the  vector  from  the  libration  point  to  the  satellite,  p,  where 

r  =  r0+p,  (2.27) 

and  p  is  expressed  in  the  local  coordinate  frame.  The  individual  coordinates  can  be 
written  as 


X  =X0+x 

(2.28) 

Y  =  Y0+y 

(2.29) 

Z  =  Z0  +  z , 

(2.30) 

11 


where  X0,  Y0,  and  Z0  are  the  coordinates  of  the  libration  point. 

So  far,  the  equations  of  motion  contain  nonlinear  terms,  all  stemming  from  the 
potential,  U.  To  simplify  for  continued  analysis,  these  equations  must  be  linearized.  To 
do  so,  substitute  Equations  2.28-2.30  into  Equations  2.24-2.26  and  perform  a  Taylor 
series  expansion  on  the  partials  of  U,  about  the  libration  point.  The  series  is  truncated 
after  the  quadratic  terms,  and  the  derivatives  are  evaluated  at  the  libration  point.  This 
gives  equations  of  motion  with  respect  to  the  libration  point  in  the  local  coordinate  frame: 


x-2(Oy -U  mX-U  XYy  =0  (2.31) 

y  +  2(ak  -  Uyy  y  -  UXYx  =  0  (2.32) 

z-U2Zz=  0,  (2.33) 


where 


and 


U  XX 

Uyy 


=  d2U 

dx2 

=  d2U 
~  BY2 


d2U 

BZ2 


r0 


U  XY 


BXBY 


(2.34) 

(2.35) 

(2.36) 

(2.37) 


The  x  and  y  motions  are  coupled,  and  the  z  motion  is  independent.  Therefore,  all  cross¬ 
term  derivatives  involving  z  disappear  and  are  not  shown  above.  Also,  for  collinear 
libration  points, 
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<3 

II 

o 

(2.38) 

The  linearized  equations  of  motion  about  a  collinear  libration  point  become 

x—lojy—U^x  —  0 

(2.39) 

.  y  +  2cai-Uny  =  0 

(2.40) 

z-U2Zz  =  0. 

(2.33) 

2.4  Lissajous  Orbits  and  Quasi-Periodic  Trajectories 

By  building  a  state  vector,  x,  consisting  of  the  positions  and  velocities,  the 

equations  of  motion  can  be  expressed  in  state-space  notation  as 

x  =  Ax, 

(2.41) 

where 

X  =  [x  y  z  X  y  zj , 

(2.42) 

x  =  [x  y  z  x  y  zj , 

(2.43) 

and 

O 

i-H 

o 

o 

o 

o' 

0  0  0  0  1 

0 

0  0  0  0  0 

l 

A  = 

U  M  0  0  0  2 (0 

o  ‘ 

(2.44) 

O 

o 

1 

V 

o 

0 

0  0  Un  0  0 

0 

There  exists  a  quasi-periodic  solution  to  the  linearized  equations  of  motion. 

The 

development  here  is  very  similar  to  that  in  Wie19.  Because  the  z  motion  is  uncoupled. 

two  characteristic  equations  exist,  the  in-plane  equation, 

A4  +(4 w-Un  -Un)X2  +UXXUYY 

=  0 

(2.45) 
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and  the  out-of-plane  equation, 


A2 -Ua  =0.  (2.46) 

The  out-of-plane  equation  is  a  simple  harmonic  oscillator.  The  eigenvalues  are 

^5,6  =  |^zz|  =  -M »  (2.47) 

where  G)z  is  the  out-of-plane  frequency.  Out-of-plane  motion  is  always  periodic.  The  in¬ 
plane  eigenvalues  are 

A IJ=±V-A+i/A2+ft2  ■  (2.48) 

and 

4.  =  .  (2-49) 

where  (Q^  is  the  in-plane  frequency, 

A  =2 a>-{U**+Un) ,  (2.50) 

2 

and 

P^-U^Uyy.  (2.51) 

Szebehely1  proves  that  U ^  is  always  positive,  whereas  U ^  and  U ^  are  always 
negative.  The  in-plane  motion  has  two  oscillatory  poles,  one  stable  pole,  and  one 
unstable  pole.  By  judiciously  choosing  initial  conditions,  the  stable  and  unstable  poles 
can  be  made  to  vanish.  Let 

Q  =  [q i  —  q6\  (2.52) 

and 

P  =  \p,  -  pj.  (2.53) 
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where  qt  is  the  normalized  right-side  eigenvector  corresponding  to  A,-  and  pt  is  the 
normalized  left-side  eigenvector  corresponding  to  A, ,  then  through  modal  decomposition 

A  =  QAPt  (2.54) 

where 


(2.55) 


and 


PTQ  =  I. 


The  solution  to  the  differential  equation  in  Equation  2.40  then  becomes 

x(f)  =  eA‘x(0)  =  '£eX'JqipJx(0) . 


Note  that  pj x(0)  is  scalar.  Setting 

Pi  x(0)  =  p\  x(0)  =  0 

leaves  only  oscillatory  poles  remaining.  These  conditions  can  be  met  if 


and 


(O 

*(0)  =— ^y(O), 

K 


y(0)  =  -ns>^x(0), 

where 

^  _  ®xy  "*■  U  xx 

20x0^ 


(2.56) 


(2.57) 


(2.58) 


(2.59) 


(2.60) 


(2.61) 
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Finally,  the  quasi-periodic  solution  to  the  linearized  equations  of  motion  for  a  collinear 
libration  point  can  be  expressed  simply  as 


x(t)  =  x(0)cos(©^f)  +  ^^sin(m^f) 

(2.62) 

y(0  =  y(0)cos(o)^t)  -»a(0)sin(G)^0 

(2.63) 

z(t)  =  z(0)cos(m  f)  +  -^^sin(n)  t) . 

t 0 , 

(2.64) 

Orbits  of  this  type  are  called  Lissajous  trajectories.  A  Lissajous  orbit  is  quasi-periodic 
because  the  ratio  of  in-plane  frequency  to  out-of-plane  frequency  is  not  an  integer. 
Barring  any  forces  other  than  the  gravity  of  the  two  large  bodies,  a  satellite  with  the 
proper  initial  conditions  in  the  circular  restricted  three-body  setting  will  follow  a 
Lissajous  orbit  forever. 
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3.  LINEAR  QUADRATIC  REGULATOR  THEORY 


3.1  Open-Loop  Dynamics 

With  no  control,  disturbances,  or  perturbations,  and  proper  initial  conditions,  a 
satellite  will  follow  the  Lissajous  orbit  described  above.  This  open-loop  trajectory  is 
useful  as  a  preliminary  reference  trajectory  for  a  spacecraft.  Generally,  the  open-loop 
dynamic  equation 

xr<f(t)  =  Axr4(t)  (3.1) 

has  the  solution 

=  <&(t-t0)xr<f(t0)  =  eA^xref  (f0),  (3.2) 

if  A  is  a  constant  matrix.  The  vector  xref  (t)  is  the  desired  reference  state  at  time  t,  and 
-/„)  is  the  open-loop  state  transition  matrix.  It  is  shown  by  Nise20  that  the  state 
transition  matrix  is  easily  approximated  by  an  infinite  series  if  A  is  constant, 

d>(t-t0)  =  eA('-'»>  •  (3.3) 

S  kl 

3.2  State  Feedback  Control 

Perturbations,  natural,  man-made,  or  deliberate,  will  require  the  satellite  be 
controlled.  With  the  added  control,  the  linear  differential  equation  becomes 

x  =  Ax  +  fiu ,  (3.4) 

where  u  is  the  control  vector  and  fi  is  the  matrix  that  maps  the  control  effort  to  the  state- 
space.  The  initial  condition  x(f0)  is  also  given.  For  this  thesis,  the  system  will  remain 
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time-invariant  (A  and  B  are  constant  matrices).  The  control  is  modeled  as  ideally  applied 
acceleration  in  the  x,  y,  and  z  directions.  Therefore, 


0  0  0 
0  0  0 
0  0  0 
1  0  0 
0  1  0 
0  0  1 


and 


M, 


w. 


(3.5) 


(3.6) 


where  u,  is  the  control  in  the  x  direction,  u2  is  the  control  in  the  y  direction,  and  u3  is  the 
control  in  the  z  direction.  The  control  vector,  u,  has  the  linear  state  feedback  form 

u  =  -K(t)x,  (3.7) 

where  K  is  a  matrix  of  gains  that  relate  the  states  to  the  control.  This  is  also  known  as  the 
control  law.  Closing  the  loop  and  substituting  back  into  the  differential  equation  yields 

x  =  (A-BK(t))x,  (3.8) 


which  has  the  solution 


x(t)  =  <S>cl(t,t0)x(t0),  (3.9) 

where  <f>cl(t,t0)  is  the  time- varying,  closed-loop  state-transition  matrix.  The  state 
transition  matrix  propagates  the  state  from  some  time  (in  this  case  t0)  to  another  time  (in 
this  case  t ).  The  state  error  is 

x(f)  =  x(f)-xr^(0.  (3.10) 
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By  differentiating  Equation  3.10,  it  can  be  shown  that  x(t)  follows  the  same  dynamic 
constraint  as  the  state.  By  superposition, 

i  =  (A-BK(t))x,  (3.11) 

and 

x(t)  =  &cl(t,t0)x(t0).  (3.12) 


Also,  the  control  law  becomes 


u  =  -K(t)x . 


(3.13) 


3.3  Continuous  Optimal  Control 
The  cost  function 

J  =  j|sr  (t)Wx(t)  +  uT  (t)Vu(t)}/t  (3.14) 

h 

is  a  scalar  number  that  represents  the  desired  performance  of  the  system.  The  first  part  of 
the  cost  function  deals  with  the  state  error.  The  matrix,  W,  weights  the  state  error  and 
must  be  positive  semi-definite.  The  second  part  of  the  cost  function  deals  with  the 
control  The  matrix,  V,  weights  the  control  and  must  be  positive  definite  and  thus 
invertible.  Altering  the  matrices  W  and  V  changes  the  interaction  between  the  control  and 
state  error.  For  this  thesis,  no  weight  will  be  given  to  the  coupling  between  states,  the 
coupling  between  controls,  or  state-control  coupling.  Hence,  W  and  V  will  remain 
diagonal. 

The  following  development  comes  from  Friedland21.  The  cost  function  can  also 
be  written  in  the  form 
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(3.15) 


/  =  J[ir(T)  U^T)^  °TS<T) 


dx 


'J_u(t)J 

The  quadratic  nature  of  the  cost  function  ensures  that  it  is  always  non-negative. 
Optimally,  both  the  state  error  and  control  would  be  Zero,  as  would  J.  Therefore,  the 
purpose  of  control  is  to  minimize  the  cost  function.  Substituting  the  control  law, 
Equation  3.13,  and  the  state-transition  matrix,  Equation  3.14,  into  the  cost  function  gives 


1  =  |(xr(r0)o;(T,(0)H'<I>d(T,(0)x(((,) 

<0 

+xr(t0)O^(T,t0)^r(T)V«(T)Oc/(T,t0)x(t0)}dT 
Extracting  the  initial  conditions  from  both  sides  of  the  integral  and  compressing, 


7  =  xr(/0)| 


J  <&Tcl  (t,  t0)ty  +  KT  (t)VA"(t)]j>c;  (t,  t0  )dz 


L  o 


X(t0)' 


Next,  define 


(3.16) 


(3.17) 


S(t0,tf)  =  +  KT  (z)VK(z))bcl(x,t0)dx ,  (3.18) 


SO 


J  =  x  (t0)S(t0,tf)x(t0). 

Since  the  initial  time,  t0 ,  is  arbitrary,  then  it  holds  for  any  time  t  that 


J  =xT  (t)S(t,tf)x(t). 


It  also  holds  that 


/  =  Jxr  (T)[W  +  KT  (x)VK{x)\(x)dx . 


(3.19) 


(3.20) 


(3.21) 


Using  the  integral-derivative  theorem 
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then 


d_ 

dx 


B(x)  B(x) 

Jp(x,t;)dt;  >=  f 

A(x)  A(x) 


dF{x&  d£ -F(x,A(x))^^-+  F(x,B(x))^^-  ,(3.22) 


dx 


dx 


dx 


—  =  -xT  (/)k  +  KT  (OVf  (0)e(0  •  (3.23) 

dt 

Also,  by  differentiating  Equation  3.20, 

—  =  kT  (t)S(t,tf  )x(0  +  xT  (t)S(t, tf  )x(0  +  xr  (i t)S (t, tf  )$(t) ,  (3.24) 

dt 

which  reduces  to 

—  =  xr  (t)[(A  -  BK(t)f  S(t,tf  )  +  S(t,tf)  +  S(t,  tf  )(A  -  BK(t))}L(t) .  (3.25) 

dt 

Setting  Equation  3.23  equal  to  Equation  3.25, 

-  (W  +  KT (t)VK(t))=  (A  -  BK(t))T  S(t,tf  )  +  S(t,tf  )  +  S(t, tf  )(A  -  BK(t)) ,  (3.26) 

or 

-S(t,tf)  =  (A-  BK(t))TS(t,  tf  )  +  S(t,  tf  )(A  -  BK(t)) +W  +  KT  (t)VK(t) .  (3.27) 
The  end  condition 

S(tf,tf)  =  0  (3.28) 

may  be  used. 

Next,  we  wish  to  find  the  gain  matrix,  K(t),  which  minimizes  the  cost  function,  J. 
Assume  that  some  optimal  gain,  K  *  (t)  .corresponds  to  the  minimum  cost  function,  J  * . 
Assume  some  other  gain,  K  *  +AK ,  corresponds  to  S  *  +AS  ,  where  AS  is  necessarily 
positive  semi-definite.  For  clarity,  all  time  indices  have  been  removed  for  the  time 
being.  By  proving  that 

J*  =  xTS*x<xT(S*  +AS)x  =  xTS*x  +  x'ASx  (3.29) 
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for  all  AS ,  then  K  *  is  truly  the  gain  that  minimizes  the  cost  function.  From  Equation 
3.27, 

-(s*+As)=(S*+AS\A- B(K  *  +A£)]+  [a  -  B{K  *  +AK)J  (5  *  +AS) 


+W  +  (K  *  +A*T)r  V(K  *  +A K) 

Subtracting  Equation  3.27  from  Equation  3.30  gives 

-AS  =  AS[A  -  B(K  *  +A/0] +[A- B(K  *  +AAT)]r  AS 

+AKtVAK  +  (K  *r  V  -  S  *  B)AK  +  AKT  (VK  *  - BTS ♦) ' 

Equation  3.31  is  of  the  form  of  Equation  3.27,  so  the  solution  is 

"<«/ >  =  fa K  +  (K*t  V -S*B)AK 

t  » 

+AKt  (YK  *  -BTS*)]0cl  (Z,t)dz 


.  (3.30) 


(3.31) 


(3.32) 


where  the  time  indices  on  S  * ,  K  * ,  and  A K  are  omitted.  Since  AS  must  be  positive 
semi-definite,  then  the  non-quadratic  terms, 

(K*tV-S*B)AK  (3.33) 


and 


must  be  zero.  Setting 


then 


A Kt(YK*-BtS*), 


K*t  V-S*B=VK*-BtS*  =  0, 


K*  =  V~'BtS*. 


Hence,  the  optimal  gain  is  computed, 

K(t)  =  V~xBTS{t,tf). 


(3.34) 


(3.35) 


(3.36) 


(3.37) 
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The  linear-quadratic-regulator  is  defined  as  the  state-feedback  control  law  of  Equation 
3.13  with  the  optimal  gain  of  Equation  3.37.  By  back  substituting  Equation  3.37  into 
Equation  3.27, 

-S(t,tf)  =  S(t,tf  )A  +  ATS(t,  tf  )  -  S(t,tf  )BV~lBTS(t,tf  )  +  W.  (3.38) 

This  equation  is  a  Riccati  equation  and  has  the  end  condition  given  by  Equation  3.28. 

3.4  The  Linear-Quadratic-Regulator 

So  far,  the  continuous  controller  depends  on  knowledge  of  a  final  time.  Often,  the 
final  time  is  unknown  or  irrelevant.  The  Riccati  equation  has  an  asymptotic  solution 
when  the  system  is  linear  time-invariant,  but  since  it  is  solved  backwards  in  time,  the 
“unchanging”  part  of  the  solution  occurs  near  t0 .  In  other  words,  S  varies  mostly  around 

tf .  Therefore,  with  an  infinite  horizon  assumption, 

limv_.5(r,r/)  =  0.  (3.39) 

In  the  limit,  the  Riccati  equation  becomes 

0  =  SA+ATS-SBV1BTS+W .  (3.40) 

This  equation  is  known  as  the  algebraic  Riccati  equation,  the  solution  to  which  is  now 
constant.  Only  the  positive  definite  solution  will  be  considered.  From  Equation  3.36,  it 
follows  that 

K=V~1BtS,  (3.41) 

where  K  is  also  constant.  The  constant-coefficient  control  law 

u(t)  =  -Kx(t)  (3.42) 

is  known  as  a  continuous  time-invariant  linear-quadratic  regulator.  Once  the  control 
gain,  K,  is  constant,  then  the  solution  to  the  closed-loop  differential  equation  becomes 
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x(0  =  <&cl(t-t0)x(t0)  =  e[A-BK^>x(t0). 


(3.43) 


3.5  Stability  of  the  Linear-Quadratic-Regulator 

Stengel22  proves  that  the  optimal  control  is  always  stable.  First,  define  a  positive 
definite  Lyapunov  function 


3(x,t)  =  xT(t)Sx(t), 


(3.44) 


where  S  is  the  positive  definite  solution  to  the  algebraic  Riccati  equation.  Stengel22 
shows  that  for  asymptotic  stability,  the  time  derivative  of  the  Lyapunov  function  must  be 


negative  definite  at  all  times.  Taking  its  derivative, 


$(x,t)=^(x,t)i(t)  =  2xr(f)S[Ax(0  +  ffu(0]. 

OX 


(3.45) 


Substituting  Equations  3.41  and  3.42  for  u,  and  reducing, 


3(x,r)  =  2xr  (t)S  [a  -  BV~1BtS  )t(0 

=  xr(07[A-RF"15rs]+  [a-BV~1Bts]  5^(0- 
=  xT  (/){SA  +  ATS  -  2SBV~lBTs}t(t) 


(3.46) 


From  the  algebraic  Riccati  equation  (Equation  3.40), 

SA  +  AT S  -2SBV~l BT S  =  -W-SBV~lBTS  ; 


(3.47) 


$(x,t)=-xT(t)ty  +  SBV~1BTsk(t). 


(3.48) 


Because  W  is  positive  semi-definite,  S  is  positive  definite,  and  V  is  positive  definite,  then 
the  time  derivative  of  the  Lyapunov  function  is  negative  for  all  times.  Hence,  the  steady- 
state  optimal  control  gain  is  stabilizing. 
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4.  DISCRETE  LINEAR  QUADRATIC  REGULATOR 


4.1  Sampling  a  Continuous  System 

Often,  a  system  is  continuous,  but  only  samples  of  the  states  and  control  are  used. 
Usually,  this  is  the  case  with  spacecraft.  Although  the  spacecraft  move  continuously,  its 
position  and  velocity  are  measured  in  discrete  intervals  and  computed  using  digital 
computers.  Some  thrusters  can  be  operated  continuously,  but  most  are  only  fired 
intermittently.  For  most  discrete  systems,  the  sampling  interval,  or  time  step,  is  held 
constant.  Now,  we  wish  to  find  some  similar  corresponding  discrete  system, 

x*+i  =  Adxk+Bduk,  (4.1) 


to  our  known  continuous  system 

x(t)  =  Ax(t)  +  Bu(t) . 


(4.2) 


The  subscript  k,  refers  to  the  vector  at  time  kT,  where  T  is  the  sampling  interval.  In  other 
words,  for  any  vector  v, 

yk+m  =  v(kT  +  mT) .  (4.3) 

Also,  the  control  vector,  u,  is  generic  and  not  necessarily  of  the  form 

u(f)  =  -Kx(t) .  (3.42) 

The  solution  to  the  differential  equation  with  generic  u  is 

t 

x(t)  =  <JH.t-t0)x(t0)+jo(t-z)Bu(r)dr.  (4.4) 

h 


Evaluating  over  one  time  interval, 

t  =  kT  +  T 


(4.5) 


and 
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o 

II 

*- 

(4.6) 

Then, 

k T+T 

x(kT  +  T)  =  ®(T)x(kT)  +  j®(kT+T-z)Bu(z)dz . 

kT 

(4.7) 

If  the  control  vector  has  the  zero  order  hold  property  (piecewise-constant),  then 

c 

II 

£ 

3 

(4.8) 

for 

kT<t<kT+T , 

(4.9) 

and  is  constant  over  the  period  of  integration.  Thus, 

~kT+T 

xM  =  <S>(T)xk  +  jo(kT  +  T-z)Bdz  uk. 

kT 

(4.10) 

Phillips  and  Nagle23  use  a  change  of  variable  to  simplify  the  integral,  let 

<7  =  kT  +  T  —  Z . 

(4.11) 

Then  the  integral  becomes 

0  T  T 

-  J  <&(o)Bdo  =  J  0(<j)Bdc  =  J  <S>(z)Bdz . 

TOO 

(4.12) 

Therefore, 

“r 

xt+1  =  ®(T)xk  +  \®(z)Bdz  uk 
.0 

(4.13) 

has  the  desired  form  of  Equation  4. 1 ,  where 

Ad  =<D(D 

(4.14) 

and 

T 

Bd  =$<f>(z)Bdz. 

0 

(4.15) 
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4.2  Discrete  Optimal  Control 

Next,  the  continuous  cost  function  is  represented  as  a  summation  of  integrals, 
rather  than  a  single  integral, 

J  =  £  T{xr  (t)Wx(t)  +  uT  (t)Vu(t)}lt ,  (4. 16) 

4=0  '* 

where  kf  corresponds  to  the  final  time,  tf  .  Assuming  piecewise-constant  control  from 

above,  then  over  the  period  of  integration 

x(t)  =  Adxk  +  Bd  uk ,  (4. 17) 


and 


u(0  =  u*. 


(4.18) 


Substituting  Equations  4.17  and  4.18  into  the  cost  function,  Equation  4.16,  and  reducing, 

J  =  1L  JfcM  +uj \BTd)W{Adxk  +  Bduk)  +  uTkVuk}it 
*=o  ,t 


y  1  kp 

=  X  \{*TkAdWAjXk  +x[AjW^u*  +uTkBTdWAdxk 
4=0  h 

WkBTdWBdnkWkVnk}dt 

X 


-‘fjfjfe  »; 
4=0  [  1 


AdWAd 


ATdWBd 


BTdWAd  (BTdWBd  +  V)j 


u, 


dt 


(4.19) 


Both  xk  and  ut  are  constant  over  the  period  of  integration  and  can  be  removed  from  the 
integral 


/-*£fe  -iff 

*=0  tL 


ATdWBd 


ATdWAd 


BldWAd  (Bd  WBd  +V\ 


\dt\ 


x. 


u. 


(4.20) 


This  can  be  rewritten  as 
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y  =  V[;[  <["'■;  ""I**], 

a  W  V,J»J 

where 

(4.21) 

fk+ 1 

(4.22) 

V,  = 

(4.23) 

and 

Md  =  JA>£>. 

(4.24) 

These  integrals  remain  constant  as  long  as  the  sampling  interval  is  constant;  therefore, 

they  can  be  expressed  as 

T 

Wd  =  j  AdWAddt , 

0 

(4.25) 

T 

Vd  =  j(BdWBd  +V)dt , 

0 

(4.26) 

and 

T 

Md  =  j  AdWBddt . 

0 

(4.27) 

Special  care  must  be  given  to  the  imbedded  integral,  Bd .  The  upper  limit  must  be 

altered  from  the  time  step,  to  the  variable  of  integration  for  the  outer  integral  when 

calculating  Vd  and  Md .  Thus, 

T  t 

Bd  =  J  <S>(t)Bdx  =  J  G>(x)Bdx . 

0  0 

(4.28) 
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The  coupled  state-control  weighting  term,  Md ,  is  interesting.  In  the  continuous 
system,  no  coupling  exists  (although  linear  quadratic  regulator  theory  could  be  modified 
to  accommodate  it).  However,  once  sampling  occurs  and  the  system  becomes  discrete,  a 
coupling  term  is  introduced. 

Next,  the  cost  function  is  modified  to  include  a  weighting  on  the  state  error  at  the 
final  time.  Also,  a  factor  of  one-half  is  added  to  simplify  the  math.  This  has  no  effect 


because  optimizing  7  would  be  the  same  as  optimizing  27. 


7  -  —xTkSkxk  +  [xj  u 

2  f  f  f  *=o  2 


riWd  Md 


k\MTd  Vd 


T5*" 

_LU*_* 


(4.29) 


Eventually,  this  final  state  error  will  go  to  zero  as  the  final  time  goes  to  infinity,  but  is 
necessary  for  the  development  of  the  theory.  Next,  a  Lagrange  multiplier  is  appended  to 
the  cost  function  inside  the  summation 


J  =  ^%Skxt/  +£4[x*r  u 


T]Wd  Md  xk 


Mi  V. „  u 


(4.30) 


+%* (Ad*k  +Bduk  -xk+1)} 

From  Equation  4.1,  the  last  term  in  Equation  4.30  is  zero,  so  the  value  of  the  cost 
function  remains  unchanged.  The  Lagrange  multiplier,  Xk ,  is  also  called  the  co-state 
vector.  The  cost  function  is  rewritten  as 


i  kr~1 

j=-iTlislili  +£k-AL, 


(4.31) 


where 


Hk  =^i 


Wd  Md 

Ml  K 


+  ^l+Md^k+Bduk), 


(4.32) 
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and  Hk  is  referred  to  as  the  Hamiltonian.  Changing  the  lower  limit  of  summation  from  0 


to  1, 


i  *,-i 

J=-iISt  xt  +H„ S,  +£K  -Alx,}.  (4.33) 

^  *=1 

Assume  a  function,  f(x1,x2,...,xn ),  exists  with  a  minimum, 

/mb  (xi  *>xi *»•  •  •*  xn  *)  •  A.  location  near  the  minimum  can  be  written  as 


f(xl,x2,...,x„)  =  fmiB(x *,x2*,...,xn*)+ 


0JC, 


&,  + 


a/ 


0X, 


Sx2  +■ 


jr 

0Jt„ 


5cn  +  (higher  order  terms) 


(4.34) 


The  first  variation  of/ is  defined  as 


5f  =  f(x1,x2,...,x„)~  fwin(x1*,x2*,...,xn*)  = 


a/ 


dtc, 


+ 


3/ 


Sx2  +• 


jr 


(4.35) 


+  (higher  order  terms) 


Assuming  the  variation  is  very  small,  or/ is  very  near  /^ ,  the  higher  order  terms  can  be 
neglected.  Reference  Bryson  and  Ho24  or  Stengel22  for  more  information  on  the  topic  of 
optimality.  Furthermore,  the  following  derivation  for  discrete  optimal  control  is  standard 
and  used  in  both  Bryson  and  Ho24  and  Stengel22.  Additionally,  this  technique  could  be 
used  to  develop  continuous  optimal  control  rather  than  the  Friedland  method  used  in 
Section  3.3. 

Setting  the  first  variation  of  the  cost  function  (Equation  4.33)  to  zero  will  provide 
necessary  conditions  for  optimality.  The  first  variation  is 
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(4.36) 


5J  = 


d-xf  S 


*k{z>kfxkf  3  A*  x* 

Z _  _ / _ f 


dxt 


dx. 


du„ 


‘fl 

k= 1 


a Hk  wkx 


dxk  dxk 


dHk 

l^+3TLAl* 
du. 


The  variation  with  respect  to  the  co-state  has  been  omitted  for  brevity.  If  it  were 
included,  it  would  simply  lead  back  to  the  original  difference  equation  constraint.  To 
ensure  that  dJ  is  zero,  all  individual  variations  within  must  be  zero.  The  initial  state 
error  vector  and  initial  control  vector  are  constant,  so  their  variations  are  inherently  zero 

8x0=  0  (4.37) 


and 


Sa0=0.  (4.38) 

Evaluating  the  partials  in  Equation  4.36  and  setting  them  to  zero  leaves  the  three  vector 
equations. 


=  xX+urMj+Af 


*+lAi  » 


(4.39) 

(4.40) 


and 


=  (4.41) 

du* 

These  equations  are  known  as  the  Euler-Lagrange  equations.  Because  the  final  time  is 
arbitrary,  it  stands  from  Equation  4.39  that 

A l=*lsk,  (4.42) 

and 
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(4.43) 


^t+i  *^*+i  x*+i » 

for  all  sampling  instances.  Substituting  Equation  4.1  into  Equation  4.43, 

K+i  =  Sk+ ,  (Adxk  +  Bd  u*  ) ,  (4.44) 

and  substituting  Equation  4.44  into  Equation  4.41, 

0  =  xTkMd  +ukVd+xTkATdSkHBd  WkBTdSk+lBd 

,  (4  45) 

=  *I(Md  +  ATdSt„Bd)  +  uTl(Vd  +  BTdS„,Bd) 

Taking  the  transpose, 


0  —  (M  d  +  Bd  Sk+1Ad  )xk  +  (Vd  +  Bd  Sk+1Bd  )ut .  (4.46) 

Note  that  W,  V,  and  Sk ,  and  in  turn  Wd  and  Vd ,  are  at  least  positive  semi-definite,  so  the 
control  can  be  found  from 

ut  =-(Vd  +  Bd  Sk+1  Bdyl(MTd  +BTdSk+1Ad)xk ,  (4.47) 

which  is  in  state-feedback  form.  The  control  law  is  then 

(4.48) 

where  the  discrete  control  gain  matrix  is 

**  =<yd  +BTdSk+lBdr(MTd  +BdSk+lAd).  (4.49) 

Substituting  the  transpose  of  Equation  4.44  into  Equation  4.40  gives 

K  =*lWd  +uTtMd  +(x[Aj  +u  TkBTd)Sk+xAd,  (4.50) 

and  transposing, 

K  =wdxk  +Mduk  +ATdSk+iAdx  +  ATdSk+lBduk.  (4.51) 

Substituting  Equation  4.42  for  the  co-state  and  combining  like  terms, 

St*k  =Wd  +ATdSk+iAd)xk  +  (Md+ATdSkHBd)uk.  (4.52) 
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Replacing  the  control  with  its  state-feedback  form  and  dropping  xk  from  all  terms  leaves 
the  discrete  Riccati  equation: 

$k  =  +  A*  Sk+lAd  —  ^ 

(MTd  +BTdSk+1Ad)T(yd+BTdSk+lBdri(MTd  +BTdSkHAd)‘ 

This  equation  is  solved  backward  in  time  with  the  known  end  condition,  Skf .  However, 

just  as  with  the  continuous  case,  S  reaches  a  steady-state  value  at  the  beginning.  The 
transient  part  of  S  occurs  near  the  final  time.  If  the  final  time  goes  to  infinity,  then  S  is 
considered  constant.  Then, 

=S,„=S.  (4.54) 

and  the  discrete  algebraic  Riccati  equation  becomes 

0  =Wd  +ATdSAd  -S-{MTd  +BTdSAd)T(Vd  +  BTdSBd)~l(MTd  +BTdSAd).  (4.55) 
Furthermore,  the  discrete  control  gain  matrix  is  now  constant 

Kd=(Vd  +BTdSBdYl(MTd  +BTdSAd ).  (4.56) 

A  continuous  optimal  control  method  has  been  found.  Once  the  system  is 
sampled,  with  a  sampling  interval  the  same  as  the  maneuver  interval,  the  continuous 
system  is  transformed  into  a  discrete  system.  Finally,  an  optimal  discrete  control  method 
has  been  built. 
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5.  FORMATION  FLYING 


5.1  Uncoupled  Satellite  Control 

So  far,  all  the  development  has  focused  on  the  control  of  only  one  satellite.  For 
many  future  missions,  multiple  satellites  will  be  required.  Adding  additional  satellites  to 
the  state-space  form  is  done  by  making  the  matrices  block  diagonal  and  appending 
additional  satellite  states  on  the  state  error  vector.  For  the  continuous  case  example  from 


Equation  4.2,  the  state  error  vector  becomes 


(5.1) 


the  control  vector  becomes 


U=tl,  U2  •••  Uj}, 


the  plant  dynamics  matrix  becomes 


(5.3) 


and  the  control-mapping  matrix  becomes 


(5.4) 


The  numerical  subscript  refers  to  the  satellite  number,  with  subscript  j  being  the  last 


satellite  included  in  the  single  condensed  equation. 
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for  all  i  and  j.  Appending  the  system  in  this  manner  allows  for  multiple  satellites,  but 
they  remain  uncoupled. 
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5.2  Coupled  Satellite  Control 

Formation  flying  requires  that  at  least  one  of  the  satellites  is  controlled  relative  to 
another  satellite.  With  the  uncoupled  system,  there  is  no  relative  control.  Wagner25 
shows  that  to  control  one  satellite  relative  to  another  the  differential  equations  of  the  two 
satellites  are  simply  subtracted.  For  example,  if  satellite  2  is  controlled  relative  to 
satellite  1,  then 

x2 -Xj  =  A(x2 -Xj)  +  5(u2 -Uj).  (5.13) 

This  is  rewritten  as 

Kei  -  Axrd  +  Bu2-Bul,  (5. 14) 

where  xrd  is  the  relative  state  vector.  By  superposition, 

Ki  =  Axrd  +Bu2  -Bu{.  (5. 15) 

If  the  state  error  vector  of  satellite  2  is  then  redefined  as  the  relative  state  error  vector 
compared  to  satellite  1,  then 

x2=Ax2  +  Bu2-Bul.  (5.16) 


Similarly  for  the  discrete  system, 

*(*+D2  —  Adxk 2  +  Bduk2  —  Bdukl .  (5.17) 

If  all  additional  satellites  are  controlled  relative  to  the  first  satellite,  then  the  continuous 
state-space  system  formed  from  Equations  4.2,  5.1, 5.2,  and  5.3  holds,  but  the  control¬ 
mapping  matrix  becomes 


B  = 


*i 

~B\  B, 


(5.19) 
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For  the  discrete  system,  Equations  4.2,  5.5,  5.6,  and  5.7  remain  the  same,  but  the  discrete 
control-mapping  matrix  changes  to 


B,= 


B 


dl 


-B„  B 


dl  di 


B. 


dl 


B 


di 


B 


di 


B 


dj 


(5.20) 


For  both  the  continuous  system  and  discrete  system,  the  additional  “drone”  satellites  are 
now  controlled  relative  to  the  first  “hub”  satellite.  The  reference  trajectory  of  the  hub 
satellite  is  whatever  stellar  trajectory  is  desired.  A  Lissajous  orbit  around  the  sun-Earth 
L2  point  is  an  example.  A  reference  trajectory  of  a  drone  satellite  is  whatever  relative 
motion  or  position  is  desired  about  the  hub  satellite.  Positioning  a  drone  satellite  at  a 
point  one  kilometer  in  the  negative  along-track  direction,  with  no  relative  velocity,  would 
have  the  desired  reference  trajectory 

xrrf2=[0  -1  0  0  0  Of.  (5.21) 
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6.  DISCRETE  KALMAN  FILTER 


6.1  Estimates  and  Measurements 

Thus  far,  the  system  assumes  perfect  mathematical  modeling  and  no  noises  or 
disturbances.  In  real  life,  however,  systems  are  not  ideal  Also,  the  true  state  (or  state 
error)  is  directly  implemented  in  the  control  scheme.  Often,  the  truth  is  unknown  and  an 
estimate  of  the  state  (or  state  error)  must  be  used.  A  Kalman  filter  blends  together  state 
estimates  and  noisy  sensor  measurements  of  the  true  states.  This  is  a  method  of  reducing 
noise  and  providing  the  best  estimate  of  the  actual  unknown  states  (or  state  errors). 

Process  noise  is  added  to  the  discrete  system 

x*+i  (61) 

where  w  is  the  random  process  noise  vector.  The  process  noise  vector  is  assumed  to  be 
Gaussian  white  noise  (zero-mean,  E[w\  =  0)  with  a  known  covariance  structure 

E[wwT\=Q,  (6.2) 

where  £[a]  is  the  expected  value  of  a  as  defined  in  Brown  and  Hwang26.  The  process 
noise  is  included  to  account  for  unmodeled  disturbances,  plant  errors,  and  plant 
nonlinearities.  The  process  noise  covariance  matrix,  Q,  is  diagonal  and  symmetric. 

The  measurements  are  expressed  as  a  linear  combination  of  the  states  (or  state 

errors): 

yt=Cxt+v,  (6.3) 

where  v  is  the  measurement  noise  vector.  If  the  satellite’s  Cartesian  position  coordinates 
are  measured  exactly,  then 
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(6.4) 


C  = 


1  0  0  0  0  0 
0  1  0  0  0  0 
0  0  1  0  0  0 


Similar  to  the  process  noise,  the  measurement  noise  is  also  assumed  to  be  Gaussian  white 
noise  with  a  known  covariance  structure 

e[wt\=R.  (6.5) 

The  measurement  noise  covariance  matrix,  R,  is  also  diagonal,  invertible,  and  symmetric. 
Measurement  noise,  or  sensor  noise,  is  included  to  account  for  the  fact  that  sensors  are 
not  ideal.  For  this  thesis,  there  is  no  correlation  between  v  and  w.  Eventually,  the 
Kalman  filter  will  be  modified  to  handle  measurements  that  are  nonlinear  functions  of  the 
state  error. 


6.2  Multiple  Satellite  Measurements 

For  multiple  satellites,  the  measurement  equation  is  appended  such  that 


y*  = 


y« 

y*2 

l»«J 


and 


(6.6) 


C  = 


A 


cj. 


(6.7) 


If  the  type  of  measurements  are  the  same  for  each  satellite,  then  C,  is  the  same  for  all  i. 


There  is  also  an  implicit  assumption  that  measurements  for  one  satellite  are  uncoupled 


39 


from  measurements  for  another.  Similar  to  the  state  errors,  however,  the  measurements 
for  the  drone  satellites  are  relative  to  the  hub  satellite.  The  hub  satellite  measurements 
are  relative  to  the  origin  of  the  coordinate  system  in  use.  If  the  circular  restricted  three- 
body  dynamics  are  in  place,  then  the  measurements  would  be  relative  to  the  libration 
point. 

6.3  Optimal  Estimate  Updating 

As  mentioned  before,  the  true  state  error  is  not  known  and  a  state  error  estimate 
must  be  used.  The  estimate  is  expressed  in  two  ways.  The  a  priori  estimate,  x* , 
represents  the  estimate  prior  to  assimilating  the  measurements  at  time  tk .  The  a 
posteriori  estimate,  xt ,  represents  the  estimate  after  assimilating  the  measurements.  The 
a  posteriori  (updated)  estimate  is  propagated  with  the  same  dynamic  model  as  the  actual 
state  error,  except  noise  is  excluded  since  it  has  zero-mean. 

**+i  =Adxt+Bdut.  (6.8) 

Propagating  the  updated  estimate  gives  the  a  priori  estimate  at  the  next  time  step. 

The  optimal  blending  of  the  measurement  with  the  a  priori  estimate  to  find  the 
updated  estimate  can  be  found  by  minimizing  another  cost  function, 

j,  -*;)+(y,  -«»)’■  R_‘(y.  -«,)],  («•«) 

where 

Pk(-)  =  E[(xk-xl)(xk-x~k)T],  (6.10) 

is  the  covariance  matrix  of  the  estimation  error  before  measurements  are  included,  and 
later, 
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pk  (+)  =  £[(xt  -  X*  )(x*  -  xt  Y  ]  , 


(6.11) 


is  the  covariance  matrix  of  the  estimation  error  after  measurements  are  included. 
Minimizing  this  cost  function  has  two  goals.  The  first  is  to  force  the  difference  between 
the  updated  estimate  and  the  a  priori  estimate  to  be  small,  meaning  the  updated  estimate 
is  close  to  the  old  estimate.  The  second  is  to  force  the  difference  between  the  actual 
measurement  and  the  expected  measurement  to  be  small.  If  both  these  conditions  are 
minimal,  then  the  estimated  state  error  would  be  very  near  the  actual  state  error.  If  the 
cost  function  is  found  to  be  zero,  then  the  estimate  is  perfect.  Taking  the  derivative  of  the 
cost  function  with  respect  to  the  updated  estimate  and  setting  it  to  zero  gives 

tr  =  (*,  -i JP,-'(-)-(y,  -CxjR-'C  =  0.  (6.12) 

OXi 

Next,  transpose  the  equation, 

0  =  p,'(-)(x,  -rt)-CTR-'( y„ (6.13) 

and  move  the  updated  estimate  term  to  the  other  side, 

(/>,-'  H  +  Ctf-'C)*,  =Pl-'(-)xt+CTR~‘  y,.  (6.14) 

Solving  for  the  updated  estimate, 

=  (/>,■'  H  +  CTR'cj'  (/>,-'  (-)x,  +  C’R'y, ),  (6. 15) 

and  adding  and  subtracting  the  term  CT R~lCx~k  from  the  equation, 

Xk  =  (p-1  (-)  +  CTR~lcY  (p-1  (-)x“  +  CTR  xCTk  -  CTR~'Cxk  +  CTR~lyk ).  (6. 16) 
This  reduces  to 

**  (6-17) 

where 
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(6.18) 


4  =  (pt_1  (-) + ctr-'c)1  CtR 

is  the  estimate  gain,  filter  gain,  or  Kalman  gain.  The  updated  estimate  is  now  expressed 
as  the  optimal  linear  combination  of  the  a  priori  estimate  and  measurements. 

6.4  The  Kalman  Gain 

If  n  is  the  number  of  states  and  m  is  the  number  of  measurements,  then  the  matrix, 
(fjT1  (-)  +  CTR~1C),  is  size  n  by  n.  Usually,  the  number  of  measurements  is  less  than  the 

number  of  states  ( m<n ).  Therefore,  it  is  computationally  easier  to  invert  a  size  m  by  m 
matrix  than  an  n  by  n  matrix.  Rather  than  just  simply  inverting  the  n  by  n  matrix,  a 
matrix  inversion  lemma  is  used  so  the  only  inversion  is  done  on  matrices  of  size  m  by  m. 
Let 

F  =  (/>  '  (-)  +  CTR~lcY ,  (6. 19) 

then, 

P~l  (-)  +  CTR~lC  =  F”1 .  (6.20) 

Pre-multiply  both  sides  by  F, 

FP~l  (-)  +  FCtR'C  =  7 ,  (6.21) 

post-multiply  both  sides  by  Pk  (-) , 

F  +  FCtR~ 1  CPk  (-)  =  Pk  (-) ,  (6.22) 

post-multiply  both  sides  byCr , 

FCt  +  FCTR~lCPk  ( -)CT  =  Pk  (■ -)CT ,  (6.23) 

and  factor  out  to  the  left  FCTR~l , 
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FCTR~l  (R  +  CPk  (-)Cr )  =  Pk  (- -)CT . 

Post-multiply  by  the  inverse  of  the  m  by  m  matrix, 

(6.24) 

FCTR~l  =  Pt  (-)CT  (r  +  CPk  (-)Cr  Y , 

and  post-multiply  by  CPk  (-) , 

(6.25) 

FCTRlCPk  (-)  =  Pk  ( -)CT  (R  +  CPk  ( -)CT  y  CPk  (-) . 

Subtracting  Equation  6.26  from  Equation  6.22  reveals 

(6.26) 

F  =  Pk  (-)  -  Pk  (~)CT  (R  +  CPk  (~)CT  y  CPk  (-) , 

or 

(6.27) 

(/>,-'  (-) + ctr-'cY  =  pt  (-){r  -  cT  (r  +  cpk  (~)CT  y  CP,  H  j. 

(6.28) 

Substituting  the  result  of  the  matrix  inversion  lemma,  Equation  6.28,  into  the 

Kalman  gain  equation,  Equation  6.18,  gives 

Lk  =  Pk  (-)} -CT{R  +  CPk  (-)Cr  CPk  . 

This  is  rewritten  as 

(6.29) 

Lk  =  Pk  (-)CT  }-(/?  +  CPk  (~)CT  y  CPk  (-)CT  . 

Factoring  out  the  inverse  term  gives 

(6.30) 

Lk=Pk  ( -)CT  (R  +  CPk  (~)CT  y  {R  +  CPk  (~)CT  )-  CPk  ir)CT  JfT‘ , 

which  reduces  to 

(6.31) 

Lk=Pk(-)CT{R  +  CPk(-)CTy. 

This  is  the  common  equation  found  in  literature  for  the  optimal  Kalman  gain. 
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(6.32) 

6.5  Estimation  Error  Covariance  Propagation 

Now  that  the  estimate  can  be  updated  and  propagated,  the  estimation  error 
covariance  matrix  must  be  updated  and  propagated  as  well.  Going  back  to  the  definition 
of  the  a  priori  estimation  error  covariance  matrix  in  Equation  6. 10,  it  stands  to  reason 
that  at  the  next  time  step 

pk+i  (“)  =  "  *kn  X**+i  “ : Y 1  •  (6-33) 


Substituting  Equation  6.1  into  Equation  6.33, 

pk+ 1 (")  =  £[(4/ x*  +  Bduk  +  w-Adxk  ~Bduk) 

(Adxk  +Bdnk  +  w-Adx*  -  Bduk)T ]’ 


which  reduces  to 

PMH  =  E[(Ad(xk -xk)  +  w)(Ad(xk  -xk)  +  w)T]. 

Expanding  yields 

Pk+1(-)  =  E[Ad(xk  -xk)(xk  -xk)TATd]  +  E[Ad(xk  -xk)wT]  + 
E[w(xk-xk)TA^]  +  E[wwT] 

Assuming  the  process  noise  and  estimation  error  have  no  correlation,  then 
E[Ad(xk  -xk)wT]  =  E[w(xk  -xk)T Ad]  =  0. 

Also,  Ad  has  no  randomness  and  can  be  removed  from  the  expectation,  so 
Pk+l(-)  =  AdE[(xk  -xk)(xk  -xk)T]Ad  +E[wr]. 
Then,  from  previous  definitions  (Equations  6.2  and  6.11), 

Pk+iH  =  AdPk(+)ATd+Q. 

This  is  the  propagation  equation  for  the  estimation  error  covariance  matrix. 


(6.34) 

(6.35) 

(6.36) 

(6.37) 

(6.38) 

(6.39) 
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6.6  Estimation  Error  Covariance  Updating 

The  definition  of  the  updated  estimation  error  covariance  matrix  is  given  in 
Equation  6. 1 1.  Substituting  Equation  6. 17  into  Equation  6. 1 1, 

/>„(+)  =  E[(x4  -x;  -My,  -«,))(*,  -rt -Lt(y„  -Cx~t))T].  (6.40) 

Next,  substitute  in  Equation  6.3  and  multiply  through, 

Pk (+)  =  E[{xk  -xk-LkCxk- Lkv - LkCxk ) 

A  T  •  (6.41) 

(xk  -±;  -LkCxk  -Lkv-LkCx~k)T] 

Combining  like  terms, 

Pk(+)  =  E[((I-LkC)(xk  -xk)-Lkv)((I-LkC)(xk  -xk)-Lkv)T],  (6.42) 

and  expanding, 

Pk(+)  =  E[(I-LkC)(xk  -xk)(xk  -i-)r(/-4C)r]- 

E[(I-LkC)(xk -xk)vTLTk]-E[Lkv(xk  -xk)T (/ - L.Cf  ]  + .  (6.43) 

E[LtwTLTk ] 

Assuming  that  the  measurement  noise  has  no  correlation  to  the  a  priori  estimation  error, 
then 

E[(I-LkC)(xt  -x~k)vTLTk]  =  E[Lkv(xk  -xk)T (/ -LkC)T]  =  0.  (6.44) 

Extracting  from  the  expectations  constant  terms  (over  the  time  step)  gives 

Pk(+)  =  C I-LkC)E[(xk -xk)(xk  -xk)T](I-LkC)T  +  LkE[wT]LTk .  (6.45) 

Using  previous  definitions  (Equations  6.10  and  6.5),  this  can  be  rewritten  as 

Pk  (+)  =  (/  -  LkC)Pk  (-)(/  -  LkCf  +  LkRLTk ,  (6.46) 

which  is  known  as  the  Joseph  formula.  Multiplying  through, 

-L,CPl(-)+LiCPl(-)CTLl  +  LlRLTl,  (6.47) 
and  substituting  the  Kalman  gain,  Equation  6.32,  gives 
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Pk  (+)  =  Pk  (-)  -  LkCPk  (-)  +  Pk  (r)CT  (R  +  CPk  (~)CT  Y  CPk  (-] )CTLTk  + 
Pk  (~)Cr  (R  +  CPk  (~)CT  Y  RlZ  - Pk  ( -)CTLTk 


(6.48) 


Factoring  Pk  (~)CT  to  the  left  and  llk  to  the  right  yields 


Pk(+)  =  PkH-LkCPkH 

+pk  (~)ct  + cpk  (~)ct  y  cpk  ( -)ct +(r+ cpk  ( -)ct  yR-i]y 

Further  factoring  inside  the  braces  reveals 
Pk(+)  =  Pk(-)-LkCPk(-)  + 

p,  (~)cT  t* + CP,  <~)CT  )r‘  (r + CP,  (-)C )-/ j \X 


.(6.49) 


(6.50) 


and 

p,  (+)  =  p,  H  -  i,CP,  (-)  +  P, (~)CT{I  -  I}L] .  (6.5 1) 

Thus, 

Pk  (+)  =  pk  H  -  hCPk  H  =  (/  -  LkC)Pk  (-)  (6.52) 

is  the  update  equation  for  the  estimation  error  covariance  matrix.  The  Joseph  Equation, 
Equation  6.46,  is  numerically  stable,  as  it  preserves  a  positive  definite  Pk  (+)  for  any 
gain.  Equation  6.52  is  correct  only  for  the  exact  optimal  gain.  Truncation  and  round-off 
errors  may  lead  to  a  non-positive  definite  Pk  (+) . 


6.7  Kalman  Filter  Algorithm 

Brown  and  Hwang26  give  the  following  algorithm  for  the  Kalman  filter: 

1 .  Enter  prior  estimate,  xk ,  and  its  error  covariance,  Pk  (-) . 

2.  Compute  Kalman  gain, 

Lk  =  Pk  (-)Cr  (R  +  CPk  ir)CT  y .  (6.32) 
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3.  Update  estimate  with  measurement, 

x*  =x*+4(yt-C£*).  (6.17) 

4.  Compute  error  covariance  for  updated  estimate, 

(+)  =  (/- Lt  C)n  (-).  (6.52) 

Here,  Equation  6.46  may  be  substituted  for  Equation  6.52  to  avoid  truncation  and  round¬ 
off  errors. 

5.  Propagate  estimate  forward, 

Kh  =  AA*  +  B<t  «*•  (6-8> 

6.  Propagate  error  covariance  forward, 

n+iH  =  4*A(+Mj+0.  (6.39) 

7.  Repeat. 

6.8  Adjusting  for  Continuous  System  Noise 

In  Equation  6.1,  noise  is  appended  to  the  discrete  system.  If  noise  is  appended  to 
the  continuous  system, 

x(0  =  Ax(f)  +  Bu(t)  +  wc  ( t ) ,  (6.53) 

where  the  Gaussian  white  noise  strength  is 

E[wc(t)wTc  (r)]  =  Qc(t)S(t- T) ,  (6.54) 

then  it  must  be  discretized  to  match  the  form  in  Equation  6.1.  From  Section  4. 1,  we 
know  the  discrete  process  noise  vector,  w,  is  determined  by 

T 

w  =  JO(r)wc(T)dT.  (6.55) 


47 


Substituting  Equation  6.55  into  Equation  6.2  gives 


T  T 

Q  =  E[wwt]  =  E  (T2 )Or (t2 )dzldz2 


Lo  o 


(6.56) 


Moving  the  expectation  inside  the  integral  and  removing  non-random  processes 


o  o 


(6.57) 


and  substituting  Equation  6.54  into  Equation  6.57  yields 

Q  =  — T2)Or(T2MT,JT2  .  (6.58) 

0  0 

According  to  Maybeck27,  this  reduces  to 

T 

Q  =  j®(z)Qc(z)<l>T(z)dz.  (6.59) 

0 

Given  the  Gaussian  white  process  noise  strength,  Equation  6.59  can  be  used  to  convert  it 
to  the  form  needed  in  the  discrete  Kalman  filter  (Equation  6.39). 


6.9  Nonlinear  Measurements 

Often,  the  exact  Cartesian  position  coordinates  cannot  be  measured  directly,  or  even 
as  a  linear  combination.  The  Kalman  filter  can  be  easily  adapted  to  account  for  such 
nonlinearities.  The  measurement  equation  changes  from  the  linear  form,  given  in 
Equation  6.3,  to  the  nonlinear  form, 

yt  =m(xk)  +  v,  (6.60) 
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where  m(xk  )  is  the  vector  containing  nonlinear  functions  of  the  state  error.  The  values 


of  m(xk )  come  directly  from  the  sensor  measurements.  An  estimated  measurement  is 
calculated  by 

y  *  =m(x“).  (6.61) 

The  estimate  update  equation  (Equation  6.17)  is  modified  to  take  into  account  the 
nonlinearities,  by 

**  =K +4(y*-y*)-  (6-62> 

Note  that  Equations  6.62  and  6. 17  are  identical  for  the  linear  case  when 

yk=Cx~k.  (6.63) 

Usually  the  state  vector  is  used  in  the  measurement  rather  than  the  state  error 

vector, 

y*=m(xt)  +  v,  (6.64) 

where 

xk  =  xk  +  xrf .  (6.65) 

The  subtraction  in  Equation  6.62  ensures  that  the  result  is  the  same.  However,  special 
consideration  must  be  given  to  the  measurement  estimate,  Equation  6.61.  The 
substitution, 

y*  =m(x*),  (6.66) 

where 

(6-67) 

must  be  made. 
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The  Kalman  gain  (Equation  6.32)  is  a  function  of  the  linear  coefficient  matrix,  C. 
In  the  nonlinear  case,  C  can  be  replaced  with  H,  where 

H=  —  (6.68) 


c/x  ■  _ 

‘  *» 


(6.69) 


H  is  the  derivative  of  the  nonlinear  vector  function,  m,  with  respect  to  the  state  or  state 
error  vector,  evaluated  at  the  a  priori  estimate.  When  H  is  calculated,  the  a  priori 
estimate  is  the  best  known  estimate  at  the  time.  For  the  linear  case. 


H=—  =  C . 

dxt  „ 

*  ** 


(6.70) 


Calculating  H  requires  taking  the  derivative  of  a  vector  with  respect  to  another 
vector.  If  there  are  two  vectors 


=  k  •"  amJ 


(6.71) 


b  =  k  -  bj. 


(6.72) 


9a, 

da. 

da, 

dbx 

db^ 

db„ 

da2 

da2 

da2 

dbx 

db2 

7 

tom 

dbx 

db2 

ton. 

(6.73) 
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Modifying  the  Kalman  gain,  Equation  6.32,  to  include  the  first-order 
approximation  of  the  nonlinearities,  results  with 

Lk  =Pk(-)HT{R  +  HPk(-)HTY.  (6.74) 

Also,  to  compute  the  updated  estimation  error  covariance  the  same  substitution  is  made 

Pk(+)  =  (I-LkH)Pk(-).  (6.75) 

6.10  Range,  Azimuth,  and  Elevation 

For  satellites,  a  ranging  system  that  gives  the  scalar  distance  between  the  satellite 
and  Earth,  or  two  satellites,  is  common.  In  addition  to  the  range,  two  angles,  azimuth  and 
elevation,  are  also  measured.  With  these  three  measurements,  the  position  can  be 
determined. 

For  the  drone  satellites,  positions  are  determined  relative  to  the  hub  satellite. 
Range  is  scalar,  so  whether  it  comes  from  a  sensor  on  the  hub  or  the  drone  is 
unimportant.  Azimuth  and  elevation,  on  the  other  hand,  usually  come  from  sensors  on 
the  drone  spacecraft.  They  are  in  the  frame  of  reference  of  the  local  coordinate  system 
centered  on  the  drone.  Because  the  states  of  the  drone  spacecraft  are  represented  with 
respect  to  the  hub  spacecraft,  and  the  angles  are  with  respect  to  the  drone  spacecraft,  a 
coordinate  transformation  is  required.  However,  if  the  local  coordinate  system  on  the 
drone  is  oriented  the  same  as  the  reference  coordinate  system  on  the  hub  (bi  is  aligned 
with  x,  bz  is  aligned  with  y,  and  b3  is  aligned  with  z),  the  position  of  the  drone  can  be 
determined  from  simple  trigonometry. 
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Figure  4.  Range,  azimuth,  and  elevation 

For  convention  positive  azimuth  is  defined  as  counter-clockwise  from  the  1)2  direction  in 
the  bi-l>2  plane.  Positive  elevation  is  defined  from  the  bi-b2  plane  upwards  in  the  positive 
b3  direction.  The  position  of  the  drone  relative  to  the  hub  is  then 

x  =  rcos(el)sin(az)  (6.76) 

y  =  -r  cos  (el)  cos (az)  (6.77) 

z  =  -rsin(ef),  (6.78) 

where  r  is  the  range,  el  is  the  elevation,  and  az  is  the  azimuth. 

For  one  spacecraft, 

m  =  [r  el  azf ,  (6.79) 

where 


r  = 


[x  y  zT 

X 

y 

=  V*2  +  y2+z2, 


(6.80) 
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(6.81) 


el  =  sin 


\ 

j 


and 


az  =  sin 


rcos(e/) 


(6.82) 


Clearly  m  is  a  nonlinear  function  of  jc,  y,  and  z.  From  Equation  6.69  and  Equation  6.73, 


H  = 


Evaluating  the  partials. 
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daz 

dx 


1 


2  2  \ 
X  z 


r  cos (el)  r  co s(el)  r  cos  (el) 


Vi 


r2  cos2  (el) 


(6.90) 


and 


daz 

dy" 


-xy 


xyz 


r3  cos(el)  r5  cos3  (el) 


1- 


r  cos (el) 


daz 

dz 


dr  dr  dr  del  del  del  daz  daz  daz 


=  0. 


(6.91) 


(6.92) 


(6.93) 


dx  dy  dz  dx  dy  dz  dx  dy  dz 
For  multiple  spacecraft,  H  is  appended  in  a  block  diagonal  fashion,  consisting  of  the  H’s 
from  the  individual  satellites. 

For  the  hub  satellite,  the  same  method  works,  assuming  that  its  local  coordinate 
frame  from  which  azimuth  and  elevation  are  measured  is  oriented  identically  to  the 
coordinate  system  the  motion  is  based  on.  For  the  circular  restricted  three-body  problem, 
the  libration  point  would  act  as  the  “hub,”  and  the  hub  spacecraft  would  act  as  a  “drone.” 
Realistically,  sensors  cannot  detect  the  libration  point,  so  the  origin  of  the  local 
coordinate  system  is  translated  in  the  negative  radial  direction  to  the  earth.  Thus,  the 
earth  acts  as  the  “hub”  and  the  hub  spacecraft  would  act  as  a  “drone.”  Because  the  local 


54 


coordinate  system  undergoes  only  translation,  the  dynamics  remain  the  same,  but  initial 
conditions  must  be  altered  appropriately. 

For  the  local  coordinate  system  centered  on  the  drone  spacecraft  to  always  be 
oriented  the  same  as  the  hub  coordinate  system,  spacecraft  attitude  is  not  considered. 
Realistically,  the  sensors  measuring  azimuth  and  elevation  would  be  fixed  at  a  certain 
location  on  the  spacecraft.  The  attitude  control  system  and  the  position  control  system 
would  be  tightly  coupled. 
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7.  STELLAR  IMAGER 


7.1  Mission  Background 

Stellar  Imager  (SI)  is  the  concept  for  a  space-based,  UV-optical  interferometer, 
proposed  by  Carpenter  and  Schrijver28, 29 .  The  purpose  of  the  mission  is  to  view  many 
stars  with  a  sparse  aperture  telescope  in  an  attempt  to  better  understand  how  stars  work. 
Hopefully,  SI  will  further  the  understanding  of  the  various  effects  of  stars’  magnetic 
fields,  the  dynamos  that  generate  them,  and  the  internal  structures  and  dynamics  of  stars. 
With  such  information,  a  model  could  be  developed  to  help  forecast  solar  activity  and 
understand  the  effects  of  stellar  magnetic  activity  on  astrobiology  and  life  in  the 
Universe.  Such  a  model  could  also  help  us  anticipate  long-term  Maunder  minima  and 
grand  maxima  occurrences  of  the  sun,  which  can  change  the  overall  global  temperatures 
causing  crop  failures.  It  could  also  help  predict  short-term  solar  activity  (flares)  that  can 
disable  satellites,  knock  out  power  grids,  increase  the  speed  of  corrosion  of  oil  pipelines, 
and  jeopardize  astronauts  from  particle  radiation. 

According  to  Rarogiewicz30,  an  interferometer  is  an  instrument  that  measures 
light  waves  by  way  of  interference  phenomena.  Light  from  a  single  source  is  “picked 
off”  by  multiple  flat  (or  spherical)  mirrors,  and  then  combined  in  such  a  way  that  they 
interfere  with  each  other.  From  this  induced  interference,  extremely  small  measurements 
are  made.  For  SI,  the  light  waves  being  utilized  are  in  the  ultra-violet  frequency  range. 

SI  has  two  primary  science  goals — imaging  of  stellar  surface  activity  and 
asteroseismology.  Currently,  the  sun  is  only  one  piece  of  data,  and  thus  provides  for 
insufficient  constraints  on  theories  of  dynamos,  turbulence,  structure,  and  internal 
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mixing.  More  future  missions  are  planned  to  find  and  image  planets  in  other  solar 
systems  using  IR-interferometry.  SI  will  study  the  central  stars  of  these  solar  systems  to 
determine  the  impact  on  the  habitability  of  the  surrounding  planets. 

7.2  Performance  Goals  and  Design  Requirements 

Carpenter29  has  outlined  three  primary  performance  goals:  “1)  Image  a 
substantial  sample  of  nearby  dwarf  and  giant  stars  representing  a  broad  range  in  magnetic 
activity,  obtaining  a  resolution  of  order  1000  total  pixels,  2)  study  a  sample  in  detail, 
revisiting  over  many  years,  and  measure  sizes,  lifetimes,  and  emergence  patterns  of 
stellar  active  regions,  surface  differential  rotation,  field  dispersal  by  convective  motions 
and  meridional  circulation,  directly  image  the  entire  convection  spectrum  on  giant  stars 
and  supergranulation,  and  3)  enable  asteroseismology,  using  low  to  intermediate  degree 
non-radial  modes  to  measure  internal  stellar  structure  and  rotation.”  The  first  goal  states 
that  the  interferometer  must  view  multiple  stars,  which  implies  that  it  can  aim  in  different 
directions.  The  second  goal  forces  the  mission  to  have  a  lifetime  of  at  least  10  years. 

For  stellar  surface  activity,  ultra-violet  images  are  necessary  to  differentiate 
between  bright  spots  (strong  magnetic  fields)  and  the  surrounding  stellar  surface. 
Integration  time  is  the  amount  of  time  the  interferometer  must  stay  focused  on  the  star  to 
avoid  smearing  due  to  rotation,  proper  motions,  and  activity  evolution.  For  surface 
imaging,  the  integration  times  will  vary  from  hours  for  dwarfs,  to  days  for  giants.  For 
asteroseismology,  the  integration  times  will  vary  from  minutes  for  dwarfs,  to  hours  for 
giants. 
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7.3  Preliminary  Mission  Design 

A  preliminary  mission  design  for  the  SI  mission  was  developed  at  the  Integrated 
Mission  Design  Center  (IMDC)  at  Goddard  Space  Flight  Center.  The  leading  concept  for 
SI  is  a  500  meter  diameter,  Fizeau  interferometer  composed  of  30  small  drone  satellites 
and  one  hub  satellite.  Each  of  these  drone  satellites  will  have  a  mirror  that  reflects  the 
incoming  light  back  to  the  hub.  The  Fizeau  type  refers  to  the  method  of  combining  the 
incoming  light  beams.  After  recombining  the  light  beams  into  useful  data,  the  hub  will 
transmit  this  information  back  to  Earth. 

7.3.1  Formation  Design 

Through  simulation,  Allen  and  Rajagopal29  determined  that  30  drones  with  a 
maximum  distance  of  500  meters  between  any  two  elements  would  provide  useful  data. 
The  following  figure  summarizes  their  efforts. 
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Figure  5.  Simulated  interferometric  images  of  a  sun-like  star  at  4  parsecs 
The  object  on  the  bottom  right  is  a  model  stellar  image  of  a  sun-like  star  4  parsecs  away. 
The  other  images  are  simulated  results  of  different  size  and  shape  interferometers 
viewing  the  model.  The  first  row  has  6  elements,  or  drone  satellites,  and  the  second  row 
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has  12.  The  satellites  are  arranged  in  a  Y-formation.  The  third  row  has  30  elements 
arranged  in  a  Golomb31  rectangle  configuration.  The  first  and  third  columns  have  a 
maximum  baseline,  or  distance  between  any  two  elements,  of  250  meters.  The  second 
and  fourth  columns  have  a  maximum  baseline  of 500  meters.  The  first  two  columns 
assume  the  formation  is  stationary  and  takes  a  snapshot  of  the  model.  The  third  and 
fourth  columns  assume  the  formation  rotates  24  times  in  15  degree  increments.  Several 
conclusions  can  be  drawn  from  the  figure.  First,  with  more  drone  satellites,  the 
interferometer  recreates  a  better  image  of  the  model  Second,  the  larger  the  maximum 
baseline,  the  better  the  model  is  recreated.  And  third,  the  more  snapshots  that  are  taken, 
based  on  rotations  of  the  formation,  the  more  accurate  the  image.  It  has  been  determined 
that  30  drone  satellites  with  a  maximum  baseline  of  500  meters  provide  enough 
resolution  for  useful  scientific  data.  However,  for  some  stars,  a  second  snapshot  may  be 
necessary,  requiring  the  formation  to  reconfigure  once  with  a  90  degree  rotation. 

A  30  element  Golomb31  rectangle  can  be  imagined  as  a  30-by-30  non-diagonal 
matrix  made  up  of  zeros  and  ones,  with  ones  representing  elements.  Every  row  and 
column  has  one  and  only  one  element  in  it.  Although  this  explanation  gives  the  image  of 
the  elements,  or  drones,  lying  in  a  plane,  they  actually  lie  on  the  surface  of  a  sphere. 

The  hub  satellite  lies  halfway  between  the  surface  of  the  sphere  containing  the 
drones  and  the  origin.  Focal  lengths  of  both  0.5  km  and  4  km  are  being  considered.  This 
would  make  the  radius  of  the  sphere  either  1  km  or  8  km.  The  following  figures  clarify 
the  two  design  options. 
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Drones 


500  m 


Sphere  Origin 


Hub  500  m 

500  m  (nottoscale) 

Figure  6.  Satellite  formation  with  focal  length  0.5  km 

Drones 


Sphere  Origin 


Hub  4  km 

500  m  (not  to  scale) 

Figure  7.  Satellite  formation  with  focal  length  4  km 

7.3.2  Orbit  Design 

The  type  of  orbit  and  location  in  space  is  an  important  part  of  mission  design.  For 
SI,  the  options  considered  were  a  low-Earth  orbit,  an  Earth-following  orbit,  and  a 
Lissajous  orbit  at  a  sun-Earth  libration  point.  The  strong  gravity  gradient  for  a  low-Earth 
orbit  causes  increased  difficulty  for  the  formation  flying  precision  necessary  for  SI.  Also, 
scattered  light  from  the  proximity  of  the  Earth  would  disrupt  the  interferometry.  Because 
of  the  length  of  the  mission  (10  years  minimum)  and  the  large  number  of  satellites  (31), 
an  individual  satellite  failure  is  probable.  Replacements  will  be  necessary  and  upgrades 
are  likely.  Re-stocking  the  formation  in  an  Earth-following  orbit  is  cost  prohibitive.  A 
similar  problem  exists  for  orbits  around  the  L3,  L4,  and  L5  sun-Earth  libration  points.  At 
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LI,  between  the  Earth  and  the  sun,  SI  would  have  to  deal  with  stray  light  coming  from 
both  the  sun  and  the  Earth.  However,  at  L2,  the  sun  and  the  Earth  are  on  the  same  side  of 
SI,  reducing  unwanted  light  and  enhancing  data  collection.  The  best  orbit  choice  for  the 
formation  is  a  Lissajous  orbit  around  the  sun-Earth  L2  point.  The  amplitude  of  the 
Lissajous  orbit  will  be  about  600,000  km,  but  is  not  critical  to  the  mission.  With  this 
orbit,  SI  will  be  able  to  cover  the  entire  sky  every  half  year  while  maintaining  an  aim 
perpendicular  to  the  sun.  For  useful  imaging,  SI  must  aim  within  10  degrees  of 
perpendicular  from  the  sun. 

7.3.3  Control  Design 

To  function  properly,  SI  will  need  to  accommodate  a  wide  range  of  control 
functions.  The  formation  must  slew  about  the  sky  requiring  movement  of  a  few 
kilometers  and  attitude  adjustments  of  up  to  180  degrees.  While  imaging,  though,  the 
drones  must  maintain  position  within  3  nanometers  of  accuracy  in  the  direction  radial 
from  the  hub  and  within  0.5  millimeters  of  accuracy  along  the  sphere  surface.  The 
accuracy  required  for  attitude  control  while  imaging  is  5  milli-arcseconds  tip  and  tilt 
(rotations  out  of  the  surface  of  the  sphere).  The  rotation  about  the  axis  radial  from  the 
hub  (rotation  within  the  sphere)  is  a  much  less  stringent  10  degrees. 

Leitner  and  Schnurr32  in  the  IMDC  propose  a  three-tiered  formation  control 
approach.  The  first  tier  is  the  rough  control  using  radio  frequency  (RF)  ranging  and 
modified  star  trackers  for  sensors  and  thrusters  for  actuators.  The  relative  positions  will 
be  controlled  to  within  a  few  centimeters.  This  level  will  drive  lost-in-space 
emergencies,  formation  initialization,  large  translations  due  to  formation  slewing, 
collision  avoidance,  and  Lissajous  orbit  maintenance.  The  modified  star  trackers  will  use 
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visual  navigation  techniques  to  initialize  and  coarsely  maintain  the  formation.  The  RF 
ranging  system  will  provide  range  measurements,  and  the  modified  star  trackers  will 
provide  azimuth  and  elevation  measurements.  The  actuators  for  this  tier  will  be  four  low- 
thrust,  high  specific  impulse  (Isp)  thrusters.  The  thrust  level  will  be  on  the  Newton  to 
milli-Newton  order  of  magnitude.  This  thesis  is  primarily  concerned  with  the  first-tier  of 
the  control  structure. 

The  second  tier  is  the  intermediate  control  with  a  modulated  laser  ranging  system 
as  sensors  and  thrusters  as  actuators.  The  relative  positions  will  be  controlled  to  within 
50  microns  at  this  tier.  This  level  will  drive  primary  attitude  adjustments  and  small 
translation  maneuvers.  Twelve  10-100  micro-Newton  Indium  Field  Emission  Electric 
Propulsion  (FEEP)  thrusters  will  be  used  for  this  level.  This  propulsion  technology  is 
currently  available,  but  by  2015  should  be  vastly  improved.  Basically,  this  tier  functions 
to  smooth  the  transition  from  the  rough  control  of  the  first-tier  to  the  fine  precision 
control  of  the  third-tier. 

The  third-tier  is  the  fine  precision  control.  At  this  level,  the  satellites  themselves 
will  not  move,  but  rather  the  mirrors  will  be  adjusted  by  extremely  accurate  mechanical 
devices  with  an  accuracy  in  the  nanometer  range.  Rather  than  having  a  traditional  sensor 
to  determine  measurements,  phase  diversity  and  wave-front  error  (WFE)  sensing 
algorithms,  using  data  from  the  incoming  light  beams,  will  determine  the  needed  control 
Currently,  phase  diversity  and  WFE  sensing  are  in  their  infancy  for  use  with  spacecraft 
and  formation  flying  concepts. 
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7.3.4  Other  Subsystem  Design  Considerations 

Current  launch  vehicles  have  the  capability  to  deploy  the  hub  and  drones  to  the 
near  L2.  The  preferred  option  from  the  IMDC  is  a  Delta  III  3940-11  to  carry  the  hub  and 
a  Delta  IV  4450-14  to  carry  all  30  drones  and  a  dispenser.  Both  the  dispenser  and  the 
hub  would  have  a  discardable  hydrazine  propulsion  system  to  correct  launch  vehicle 
errors  and  aid  in  insertion  into  the  Lissajous  orbit.  Initial  total  mass  assumptions  are  100 
kg  for  each  drone  satellite,  550  kg  for  the  hub,  and  480  kg  for  the  dispenser.  The  launch 
vehicles  mentioned  can  easily  accommodate  such  masses  with  plenty  of  margin. 

The  hub-drone  setup  lends  itself  to  a  centralized  method  of  communications  and 
data  handling.  The  drones  would  talk  to  the  hub,  and  the  hub  would  process  information 
and  talk  to  the  Earth.  Relaying  the  information  from  the  hub  to  the  Earth  would  be 
enhanced  by  a  central  communications  satellite  at  L2  for  all  missions  flying  there.  The 
computer  processing  power  required  for  a  centralized  approach  would  be  immense,  but 
may  not  be  a  problem  by  2015.  Estimated  power  requirements  could  be  handled  by 
body-mounted  solar  arrays  on  both  the  drones  and  the  hub.  The  only  challenge  for  the 
thermal  system  would  be  holding  the  mirror  temperatures  constant.  With  precision 
needed  at  the  nanometer  level,  it  is  imperative  to  keep  the  mirrors  isothermal  to  avoid 
contraction  or  expansion  due  to  heating. 
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8.  SIMULATION 


8.1  Background 

An  essential  design  consideration  for  Stellar  Imager  is  how  much  on-board  fuel  is 
required  to  complete  the  mission.  This  thesis  focuses  on  the  preliminary  analysis  of  the 
position  control  of  the  formation  to  determine  the  fuel  requirements.  All  spacecraft  are 
considered  “black  boxes,”  and  satellite  attitudes  are  neglected.  All  satellites  in  the 
formation  are  considered  identical,  so  Equations  5.9-5.12  are  valid.  Also,  the  drones  are 
all  controlled  with  respect  to  the  hub,  so  the  coupled  satellite  control  equations  in  Section 

5.2  are  used.  Three  different  scenarios  make  up  the  position  control  problem — 
maintaining  the  Lissajous  orbit,  slewing  the  formation  to  aim  at  another  star,  and 
reconfiguring  the  formation  to  take  another  snapshot  of  a  star  when  necessary.  These 
three  scenarios  are  treated  independently. 

8.1.1  Earth-Sun  L2  Circular  Restricted  Three-Body  Dynamics 

Specific  numerical  values  are  needed  to  build  the  dynamics  matrix  in  Equation 
2.44  for  motion  of  a  satellite  in  the  vicinity  of  the  sun-earth  L2  point.  Table  1 
summarizes  the  constant  parameters  in  the  sun-earth  system,  given  by  the  Astronomical 
Almanac33. 
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Parameter 

Equation 

Value 

Mass  of  the  sun  (msun) 

N/A 

1.9891e30  kg 

Mass  of  the  earth  (mearth) 

™sun 

332946.0 

5.9742e24  kg 

Mass  of  the  moon  (mmoon  ) 

meanh  *0.01230002 

7.3483e22  kg 

Universal  Gravitation  Constant 
(G) 

N/A 

0.4980621312 

m 

kg -day2 

Astronomical  Unit  (AU  and  D ) 

N/A 

1.49597870ell m 

Distance  from  system 
barycenter  to  sun  (Dj ) 

m  ,  4-  m 

*earth  moon  *  \\J 

msun  +  m'anh  +  mrnoon 

4.54841e5  m 

Distance  from  system  barycenter 
to  the  earth-moon  barycenter 

(02) 

ms«n  *  AU 

msun  +  m earth  +  mmoon 

1.49597ell m 

Angular  velocity  of  system  (coi) 

lG*(msun  +  mearth+mnwon) 

V  AU 3 

(Equation  2.6) 

rad 

0.017202  — 
day 

Table  1.  Constant  parameters  of  the  sun-earth  circular  system 


Based  on  the  angular  velocity  of  the  system,  the  period  of  rotation  comes  out  to  be  365.26 
days,  or  one  year,  as  expected. 

Szebehely1  gives  the  location  of  the  L2  point  (which  he  calls  LI)  for  the  sun-earth 
system.  Note  that  the  V  direction  of  his  rotating  coordinate  frame  is  opposite  of  mine, 
leading  to  the  negative  value  of  X0 .  From  this  information,  the  partial  derivatives  in 
Equations  2.34-2.36  can  be  calculated.  The  results  are  summarized  in  Table  2. 
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Parameter 

Value 

AU* -1.0100701938  m 

Yo 

0 

Zo 

0 

U*x 

0.00263065454908155  rad. 

day 

Uyy 

-0.000871456562957885  rad2 

day 

Un 

-0.001 16737037068138  rad \ 

day 

Table  2.  Constant  parameters  of  the  sun-earth  L2  point 
Using  the  above  values,  the  dynamics  matrix  can  now  be  determined 


0  0 
0  0 
0  0 
Un  0 

0  u„ 
0  0 


0  10  0 

0  0  10 

0  0  0  1 

0  0  2co  0 

0  -Txo  0  0 

Un  0  0  0 


(2.44) 


Additional  parameters  are  needed  to  determine  the  quasi-periodic  reference  orbit 
described  by  Equations  2.62-2.64.  Table  3  lists  these  parameters. 
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Parameter 

Value 

In-plane  frequency  (o)^  ) 

0.0354038676524272  rad 

day 

Out-of-plane  frequency  (cdz) 

0.0341668021723043  rad 

day 

Non-oscillatory  poles  nulling  factor  (k) 

3.18878901709247 

Table  3.  Parameters  for  quasi-periodic  orbit  about  sun-earth  L2 


The  in-plane  period  is  177.47  days,  and  the  out-of-plane  period  is  183.90  days.  Every 
year,  the  formation  makes  roughly  two  revolutions  both  in-plane  and  out-of-plane. 

The  reference  initial  conditions  listed  in  Table  4  correspond  to  the  quasi-periodic 
Lissajous  reference  trajectory  for  one  year  shown  in  Figure  8. 


Initial  condition 

Value 

M°) 

0km 

>v(0) 

600000  km 

**( 0) 

0  km 

^(0) 

6661.563 

day 

0) 

o  *"• 

day 

-20500.081  — 
day 

Table  4.  Reference  initial  conditions  corresponding  to  Figure  8 


67 


Reference  Orbit 


x  10 


Figure  8.  One  year  quasi-periodic  reference  orbit  around  L2 


In  Figure  8,  the  origin  is  the  L2  point,  the  red  circle  is  the  starting  reference  position,  and 
the  red  ‘x’  is  the  ending  reference  position  after  365  days.  The  scale  is  in  kilometers. 

8.1.2  Generator 

So  far  the  dynamics  of  a  satellite  in  the  vicinity  of  the  sun-earth  L2  point  have 
been  approximated  with  the  circular  restricted  three-body  assumptions.  These 
assumptions  only  account  for  gravitational  forces  from  the  sun  and  Earth.  The  moon  is 
also  included,  but  not  as  an  independent  body.  The  masses  of  the  earth  and  moon  are 
combined  and  assumed  to  be  at  the  earth-moon  barycenter.  The  motion  of  the  sun  and 
the  earth-moon  barycenter  is  also  assumed  to  be  circular  around  the  system  barycenter. 

Realistically,  the  earth  has  a  slightly  elliptic  orbit  around  the  sun.  The  effect  of 
the  moon’s  gravity  will  vary  depending  on  the  exact  position  of  the  moon.  Gravitational 
forces  from  other  planets  in  the  solar  system  affect  the  motion  of  satellites  as  well.  The 
most  significant  non-gravitational  force  at  the  sun-earth  L2  point  is  solar  radiation 


pressure,  ignored  in  the  circular  restricted  three-body  problem.  Dr.  Kathleen  Howell34, 
from  Purdue  University,  has  a  program  named  Generator  that  creates  much  more  realistic 
Lissajous  orbits  than  those  derived  from  the  circular  restricted  three-body  problem. 

Using  ephemeris  files,  Generator  takes  into  account  the  effects  of  eccentricity,  an 
independent  moon,  the  other  planets  of  the  solar  system,  and  solar  radiation  pressure. 

The  resulting  Lissajous  orbit  can  then  be  used  as  a  more  accurate  reference  orbit. 

Another  benefit  of  using  a  Generator  derived  reference  orbit  is  that  the  origin  of  the 
rotating  reference  frame  is  based  on  the  earth  instead  of  the  libration  point.  Thus,  sensors 
can  be  modeled  to  measure  range  and  angles  from  the  hub  satellite  to  the  earth,  rather 
than  from  the  hub  to  the  imaginary  libration  point. 

With  the  initial  conditions  listed  in  Table  5,  the  Generator  reference  orbit  is 
shown  in  Figure  9. 


Initial  condition 

Value 

M°> 

1457251.466  km 

3W«>) 

572731.249  km 

M°) 

71916.835  km 

6779.986  *"* 
day 

S’ ref  (0) 

-5701.628  *”* 
day 

Zrefi 0) 

79.552  — 
day 

Table  5.  Reference  initial  conditions  corresponding  to  Figure  9 
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x  -Terence  for  hub  satellite 


x  10s 


Figure  9.  Generator  based  359  day  reference  orbit  around  L2 


The  origin  in  Figure  9  is  the  earth.  The  reference  orbit  begins  at  the  red  circle  and,  after 
359  days,  ends  at  the  red  ‘x\  The  scale  is  in  kilometers.  The  simulation  length  of  359 
days  is  used  to  synchronize  the  simulation  time  with  the  output  of  Generator. 

While  using  a  Generator  based  reference  orbit,  the  dynamics  of  Equation  2.44 
(based  on  the  circular  restricted  three-body  problem)  do  not  hold.  However,  in  addition 
to  providing  the  reference  positions  and  velocities,  Generator  also  numerically  computes 
and  outputs  the  dynamics  matrix,  A,  for  a  single  satellite  at  each  epoch.  Anderson 
explains  the  transformation  of  the  dynamics  matrix  from  the  inertial  reference  frame  used 
in  Generator  to  the  outputted  Earth-centered  rotating  frame  used  in  the  simulation. 

Special  consideration  must  be  given  to  the  units  of  length  and  time.  Generator  uses 
velocities  with  units  of  meters  per  second.  These  must  be  converted  to  units  of 
kilometers  per  day  for  use  in  the  simulation.  Also,  the  time  unit  within  the  dynamics 
matrix  must  be  converted  from  seconds  to  days. 
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8.2  Lissajous  Orbit  Simulation 

Following  the  Lissajous  orbit  is  not  a  problem  of  formation  control,  but  rather  a 
problem  of  maintaining  an  orbit.  Therefore,  only  the  hub  satellite  needs  simulation  to 
determine  the  amount  of  control  and  fuel  needed  to  maintain  a  Lissajous  orbit.  The 
results  can  be  extended  similarly  to  other  satellites  in  the  formation.  The  simulation  does 
permit  altering  the  number  of  satellites,  but  for  each  satellite,  the  time  it  takes  to  run  the 
simulation  increases  proportionally. 

8.2.1  Lissajous  Orbit  Simulation  Development 

First,  the  reference  orbit  is  loaded  and  converted  to  the  proper  units  from  the 
Generator  output.  The  Matlab  file  used  to  change  the  Generator  reference  orbit  output  to 
a  form  used  in  Matlab  is  given  in  Appendix  A  Next,  the  initial  conditions  are  defined. 
For  this  research,  the  satellite  is  assumed  to  start  with  no  position  or  velocity  errors.  This 
means  that  the  satellite  initial  conditions  are  identical  to  the  reference  initial  conditions 
given  in  Table  5.  In  the  simulation,  the  time  step  or  maneuver  interval,  T,  can  varied  by 
integer  values  of  days,  but  is  set  at  one  day  here. 

At  the  beginning  of  each  epoch,  the  corresponding  single  satellite  dynamics 
matrix  is  taken  from  the  Generator  output  and  converted  to  proper  units.  The  Matlab  file 
used  to  change  the  Generator  dynamics  output  to  a  form  used  in  Matlab  is  given  in 
Appendix  B.  The  continuous  control-mapping  matrix  used  is  the  same  as  Equation  3.5. 
The  state  transition  matrix  is  approximated  by  Equation  3.3,  truncating  after  the  quadratic 
term.  The  continuous  state  weighting  matrix  is  chosen  to  be 
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Ie6 


(8.1) 


and  the  continuous  control  weighting  matrix  is  chosen  to  be 


V*  1 


The  strength  of  the  process  noise  is  chosen  as 


le—6 


le  —  6 


le-6 


The  strength  of  the  process  noise  is  set  at  a  value  large  enough  to  be  noticed,  but  not  so 
much  as  to  constrict  or  destabilize  the  system.  The  system  is  then  discretized  using 
Equations  4.14,  4.25, 4.26,  4.27,  and  4.28.  The  discrete  process  noise  covariance  is 
calculated  by  Equation  6.59.  The  measurement  noise  covariance  matrix  for  the  hub 
satellite  is  chosen  to  be 


1500000 


1500000 
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The  first  term  in  the  measurement  noise  covariance  matrix  assumes  range  measurements 
from  the  earth  to  the  hub  within  0. 1  km  or  100  meters.  The  second  and  third  terms 
assume  that  the  arc  lengths  corresponding  to  the  azimuth  and  elevation  angles  are  three 
times  less  accurate  than  the  range  measurements.  To  find  the  angular  accuracy,  simply 
divide  the  accuracy  of  arc  length  by  the  range.  The  L2  point  is  approximately  1.5  million 
kilometers  from  the  earth  and  is  used  as  a  rough  constant  range. 

If  multiple  satellites  are  simulated,  then  Equations  5.7, 5.20,  and  6.7  are 
implemented.  Also  for  multiple  satellites,  the  discrete  state  weighting  matrix,  discrete 
control  weighting  matrix,  discrete  state-control  couple  weighting  matrix,  and  discrete 
process  noise  covariance  matrix  are  appended  in  block  diagonal  fashion— one  block  for 
each  satellite,  including  the  hub.  The  measurement  noise  covariance  matrix  is  also 
appended  block  diagonally,  but  the  covariance  for  the  drone  satellites  is  different  from 
that  of  the  hub.  For  drone  satellites, 

0.0003  Y 
0.5  , 

Here,  the  range  measurement  is  assumed  to  be  accurate  to  within  0.1  meters,  with  three 
times  less  accurate  arc  lengths.  The  range  from  the  hub  to  the  drones  is  the  focal  length 
of  the  interferometer  (either  0.5  or  4  km). 

The  control  gain  is  calculated  using  the  ‘dlqr’  command  in  Matlab,  which  solves 
Equations  4.55  and  4.56.  The  control  is  calculated  by  Equation  4.48,  and  the  state  error 
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and  error  estimate  are  propagated  by  Equations  6.1  and  6.8  respectively.  The  initial 
condition  on  the  error  estimate  is  set  equal  to  the  initial  condition  on  the  state  error. 

Equations  6.64-6.67  are  used  to  process  the  measurements  and  measurement 
estimates.  The  measurements  are  assumed  to  be  range,  azimuth,  and  elevation,  so 
Equations  6.79-6.82  are  valid.  The  “true”  measurements  would  realistically  come  from 
sensors,  but  for  the  simulation  purposes,  they  are  the  known  range,  azimuth,  and 
elevations  with  appended  noise.  Next,  Equation  6.83  (filled  in  by  Equations  6.84-6.93)  is 
calculated.  Finally,  the  Kalman  filter  algorithm  is  implemented  with  Equations  6.39, 
6.74,  6.62,  and  6.75  respectively.  The  initial  estimation  error  covariance  for  the  hub 
satellite  is 


1 


/»,(+)  = 


86.42 

86.42 

86.42 


(8.6) 


This  assumes  about  1  km  accuracy  for  position  (10  times  greater  than  measurement  noise 
covariance)  and  1  meter  per  second  for  velocity.  The  86.4  term  is  the  conversion  to 
kilometers  per  day  from  meters  per  second.  For  drone  satellites,  the  tolerance  is  much 
more  accurate  due  to  the  proximity  from  the  drones  to  the  hub  as  compared  to  the  hub  to 
the  earth.  Assuming  position  accuracy  of  1  m  and  velocity  accuracy  of  1  millimeter  per 
second,  the  initial  estimation  error  covariance  for  drone  satellites  is 
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.0012 

.0012 

.0012 

.08642 

.08642 

.08642 

Finally,  the  whole  process  is  repeated  at  each  epoch,  starting  with  loading  the  appropriate 
continuous  dynamics  matrix  from  Generator.  The  Lissajous  orbit  Matlab  simulation  is 
given  in  Appendix  C. 

8.2.2  Lissajous  Orbit  Simulation  Results 

The  simulation  of  a  satellite  maintaining  a  Generator  based  Lissajous  orbit  around 
L2  provides  three  main  results.  The  first  is  how  well  the  satellite  stayed  its  course,  the 
second  is  how  well  the  Kalman  filter  performed,  and  the  third  is  how  much  fuel  was 


necessary. 

Figure  10  shows  the  tracking  errors  plotted  against  time  for  one  simulation. 


hub  tracking  error 


time  (days)  time  (days) 

Figure  10.  Hub  tracking  errors 
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The  plots  on  the  left  are  position  tracking  errors,  and  the  plots  on  the  right  are  velocity 
tracking  errors.  Running  many  simulations,  it  can  be  seen  that  the  steady-state  position 
tracking  errors  are  within  0.25  km  for  each  direction,  and  the  steady-state  velocity 
tracking  errors  are  within  7.5e-4  meters  per  second  for  each  direction. 

Figure  1 1  shows  the  estimation  errors  of  the  hub  satellite  for  one  simulation. 


time  (days)  time  (days) 

Figure  11.  Hub  estimation  errors  (one  simulation) 

The  plots  on  the  left  are  the  position  estimation  errors,  and  the  plots  on  the  right  are  the 
velocity  estimation  errors.  The  blue  lines  represent  the  actual  estimation  errors,  and  the 
red  lines  represent  the  three-sigma  value  of  the  covariance.  For  any  epoch,  k,  the  three- 
sigma  value  of  the  covariance  is  calculated  by 


3ob=±3jPt(+M) 

(8.8) 

3^=±3Vn(+)(2,2) 

(8.9) 

3CTfe=±3A/Ft(+)(3,3) 

(8.10) 

3au  =  ±3-J  Pt  (+)(4,4) 

(8.11) 
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3<t#  =±3Vn(+)(5,5) 
3<tk  =±3VA(+)(6,6) . 


(8.12) 


(8.13) 


Because  the  estimation  errors  are  based  on  randomness,  many  simulations  must  be  run  to 
determine  useful  results.  Figure  12  shows  the  estimation  errors  of  a  dozen  simulations. 


Hub  Estimation  Errors  (12  simulations) 


x  lo"1 


Figure  12.  Hub  estimation  errors  (12  simulations) 

The  estimation  errors  lie  within  the  three-sigma  values  of  the  covariance  with  very  few 
exceptions.  The  position  estimation  errors  are  within  250  meters,  and  the  velocity 
estimation  errors  are  within  2e-4  meters  per  second. 

Figure  13  shows  the  control  effort  over  the  course  of  one  simulation. 
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Control  1 


0  50  100  150  200  250  300  350  400 

time 


Figure  13.  Hub  control  effort 

Remember  that  the  control  is  the  applied  directional  accelerations.  To  determine  the 
amount  of  fuel,  the  velocity  change,  or  AV ,  is  needed.  The  AV  in  each  direction  is 
found  by  numerical  integration: 


AV,=XklT 

(8.14) 

359  , 

k 

n 

< 

(8.15) 

359 

AV.=£klr> 

*=1 

(8.16) 

where  T  is  the  maneuver  interval.  The  absolute  value  of  the  control  is  taken  because  the 
direction  of  the  maneuver  has  no  bearing  on  the  fuel  used.  The  total  AV  for  one 
simulation  is  calculated  by 

AV  =  AVX  +  AVy  +  AVZ .  (8.17) 
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Averaging  the  determined  AV  from  a  dozen  simulations,  the  AV  required  to  keep  a 
satellite  in  a  Lissajous  orbit  about  L2  for  359  days  is  approximately  0.38  meters  per 
second.  Humble36  gives  the  equation  relating  propellant  mass  to  AV  as 


™prop  =  ™0 


(  z*L  ^ 


(8.18) 


V  J 

where  mprop  is  the  propellant  mass,  m0  is  the  initial  spacecraft  mass,  Isp  is  the  specific 

impulse  of  the  thruster,  and  g  is  the  gravitational  acceleration  constant  9.81  meters  per 
second  squared.  From  the  Stellar  Imager  IMDC,  Asato37  gives  the  initial  spacecraft  mass 
to  be  550  kg  for  the  hub  and  100  kg  for  a  drone.  Also,  the  Isp  for  the  low-thrust,  high-Isp 
thruster  is  assumed  to  be  10000  seconds.  Therefore,  the  amount  of  propellant  needed  to 
maintain  a  Lissajous  orbit  for  359  days  is  less  than  2.2  grams  for  the  hub  and  less  than  0.4 
grams  for  each  drone. 


8.3  Formation  Slewing 

A  key  part  of  the  SI  mission  is  to  image  many  stars.  Following  a  Lissajous  orbit 
around  L2,  SI  could  view  the  entire  sky  approximately  every  half-year  while  slewing 
about  just  the  radial  (x)  axis.  This  will  also  maintain  the  aiming  angle  perpendicular  to 
the  sun.  The  formation  slewing  simulation  follows  a  similar  algorithm  as  the  Lissajous 
orbit  simulation. 

8.3.1  Formation  Slewing  Simulation  Development 

First,  the  number  of  satellites  to  be  simulated,  the  maneuver  interval  in  minutes, 
and  the  simulation  length  in  hours  are  specified.  Also,  variable  are  the  place  in  the 
Lissajous  orbit,  ranging  from  the  first  day  to  the  358th,  and  the  slew  angle.  For 
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continuity,  the  maneuver  interval  will  be  set  to  1  minute  and  the  length  of  the  simulation 
to  24  hours.  These  24  hours  will  occur  from  the  second  to  the  third  day  as  modeled  by 
Generator.  Different  slew  angles  will  be  examined. 

The  reference  Lissajous  orbit  and  dynamics  are  loaded  from  Generator,  based 
upon  the  desired  location  in  the  orbit,  and  converted  to  the  proper  units.  Similar  to  the 
Lissajous  orbit  simulation,  the  continuous  control-mapping  matrix  used  is  Equation  3.5. 
The  state  transition  matrix  is  approximated  by  Equation  3.3,  truncating  after  the  quadratic 
term.  The  continuous  state  weighting  matrix  is  given  in  Equation  8.1,  and  the  continuous 
control  weighting  matrix  is  given  in  Equation  8.2.  The  strength  of  the  process  noise  is 
different  from  the  Lissajous  orbit  simulation, 


0 


0 


Qc  = 


o 

le-24 


(8.19) 


le-24 

le-24 

The  strength  of  the  process  noise  is  set  at  a  value  large  enough  to  be  noticed,  but  not  so 
much  as  to  constrict  or  destabilize  the  system.  The  system  is  then  discretized  using 
Equations  4.14,  4.25, 4.26, 4.27,  and  4.28.  One  key  difference  between  the  formation 
slewing  simulation  and  the  Lissajous  orbit  simulation  is  the  maneuver  interval  or  time 
step.  In  the  Lissajous  orbit  simulation,  the  time  step,  T,  is  in  days.  In  the  formation 
slewing  simulation,  the  time  step,  dT,  is  in  minutes,  and  must  be  converted  to  days  by 
dividing  by  1440.  This  is  then  used  for  determining  the  state  transition  matrix  and 
evaluating  Equations  4.14,  4.25, 4.26,  4.27,  and  4.28.  The  discrete  process  noise 
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covariance  is  calculated  by  Equation  6.59.  The  measurement  noise  covariance  matrix  is 
given  by  Equation  8.4  for  the  hub  and  Equation  8.5  for  the  drones. 

The  initial  conditions  of  the  Stellar  Imager  formation  are  coded  into  the 
simulation.  The  assumption  is  that  the  satellites  start  perfectly  in  formation,  and  slew  to 
the  new  formation  over  the  course  of  a  day.  The  new  slewed  formation  is  the  desired 

reference.  The  geometry  of  both  the  initial  formation  and  slewed  formation  are  explained 
in  detail. 

To  build  the  initial  formation,  the  hub  spacecraft  is  assumed  to  be  at  the  center  of 
a  sphere,  on  which  the  drones  lie.  The  x  direction  of  the  hub-centered  Cartesian 
coordinate  system  is  parallel  to  and  in  the  same  direction  as  the  radial  axis  of  the  earth- 
centered  rotating  coordinate  system  used  by  Generator.  The  y  direction  of  the  hub- 
centered  system  is  parallel  to  and  in  the  same  direction  as  the  along-track  direction  of  the 
earth-centered  system.  The  z  direction  of  the  hub-centered  system  is  parallel  to  and  in  the 
same  direction  as  the  cross-track  direction  of  the  earth-centered  system.  Figure  14  shows 
the  relationship  between  the  two  coordinate  systems. 
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From  Figures  6  and  7,  we  see  that  the  hub  is  not  actually  at  the  center  of  the  sphere,  but 
that  will  be  taken  into  account  later.  The  drones’  positions  are  expressed  in  the  spherical 
coordinates,  r,  6 ,  and  0 ,  which  relate  to  the  hub-centered  Cartesian  coordinates  by 

x  =  rsin(0),  (8.20) 

y  =  rcos(0)sin(0),  (8.21) 

and 

z  =  rcos(0)cos(0).  (8.22) 

The  radial  coordinate,  r,  is  the  distance  from  the  hub  to  the  drone  and  always  equals  the 
radius  of  the  sphere,  which  is  twice  the  desired  focal  length.  The  0  coordinate  is 
measured  from  the  positive  z  axis  toward  the  positive  y  axis.  Standing  on  the  positive  x 
axis  and  looking  back,  a  positive  rotation  is  clockwise.  The  6  coordinate  is  measured 
from  the  position  in  the  y-z  plane  along  the  sphere  in  a  clockwise  direction  when 
considered  from  the  positive  z  axis  looking  back. 

Initially,  the  center  of  the  formation  will  be  placed  directly  behind  the  hub  in  the 
along-track  direction.  Thus, 

#1,  =0  (8.23) 

and 

=■ y,  (8-24) 

where  6°mt  and  0  1,  are  initial  coordinates  corresponding  to  an  imaginary  central  drone 
satellite.  There  are  30  drones  for  the  SI  mission,  so  15  have  a  0°  coordinate  greater  than 
0l, »  and  15  have  a  0°  coordinate  less  than  0C°„, .  Knowing  that  the  maximum  distance 
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between  drones  is  0.5  km,  or  0.25  km  to  the  imaginary  central  drone,  the  coordinates  can 
be  calculated  by 


.0,0  .  17-1/0.25^ 
»«  =9  cm,  —  sm  - 

15  r  j 


(8.25) 


where  i  represents  the  2nd  through  16th  satellite  (the  hub  being  number  1),  and 


,o  ,o  i—16  .  -i 

0,  -^cm, - —  SU1 


f  0.25 


v  r  j 


(8.26) 


where  i  represents  the  17th  through  31st  satellite.  The  6°  coordinates  can  be  determined 
in  a  similar  fashion,  for  i  from  2  to  31, 

i  (  0.25  ) 


e,°=C,+^ta- 


v  r  ) 


(8.27) 


where  j  is  an  integer  between  -15  and  15,  such  that  each  satellite,  i,  has  a  unique  j.  This 
formation  is  a  Golomb31  rectangle  laid  out  on  a  spherical  surface  rather  than  a  plane.  The 
spherical  coordinates  for  each  drone  and  the  imaginary  central  drone  are  intermediately 
transformed  to  Cartesian  coordinates  by  Equations  8.20-8.22,  yielding  mxf,  myf ,  mzf , 

xLu »  ylem  *  and  Zce„,  •  Finally,  the  satellites  are  translated  by  half  of  the  imaginary  central 
coordinates  to  account  for  the  hub  being  halfway  between  the  origin  of  the  sphere  and  the 
surface  of  the  sphere: 


_ int^.0  _ cent 


(8.28) 


0 _ int  0  _  y  cent 

y  >  ~  y  i  „ 


(8.29) 


_0 _ int  _0  Zcent 

Zi  -  zi  — 


(8.30) 


The  drones  initially  start  with  zero  relative  velocity  to  the  hub,  so  for  i  from  2  to  31, 
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•  0  -0  •  0  n 

*,  -yt  -zt  =o. 


(8.31) 


The  reference  states  are  built  by  rotating  the  formation  through  a  slew  angle,  a . 
The  rotation  is  about  the  hub-centered  x  axis.  Thus,  for  i  from  2  to  31, 


(8.32) 


ef  =0,°, 


(8.33) 


<br *  —(b0  —cl  , 
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(8.34) 


01=01,=  0. 


(8.35) 


These  spherical  reference  coordinates  are  plugged  into  Equations  8.20-8.22,  which  yields 


the  intermediate  inlxf ,  yf ,  inl  zf ,  xrfmt ,  yTfmt ,  and  z't,  •  Just  as  in  Equations  8.28- 


8.30,  the  drones  must  be  translated  by  half  of  the  imaginary  central  drone  position: 


ref _ int  ref  _  xcent 

A:  A: 

‘  '  2 


(8.36) 
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(8.38) 


No  relative  velocities  between  the  hub  and  drones  are  desired,  so  for  i  from  2  to  31, 


x?  =y?  =z?  =  0. 


(8.39) 


The  drones’  reference  states  do  not  change  over  the  course  of  the  simulation. 

The  hub  satellite’s  initial  and  reference  conditions  are  treated  differently  than  the 
drones’,  because  the  hub’s  reference  states  change  over  the  duration  of  the  simulation. 


From  Generator,  the  hub  reference  conditions  are  known  at  the  beginning  of  the 
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simulation,  x,r<f  0 ,  and  at  the  end  of  simulation,  one  day  later,  x[efend .  The  initial 
conditions  are  set  equal  to  the  known  starting  reference  conditions, 

x?  =  x*°.  (8.40) 

The  hub  reference  states  at  some  epoch,  k,  during  the  simulation  are  determined  by  linear 
interpolation, 

xf  =  {x[efend  - xf°  +  xr*° .  (8.41) 

The  control  and  estimation  algorithm  from  one  epoch  to  another  are  identical  to 
the  Lissajous  orbit  simulation  detailed  at  the  end  of  Section  8.2.1.  The  lone  difference  is 
that  the  dynamics  are  assumed  to  be  time-invariant  over  the  course  of  the  simulation 
rather  than  being  updated  at  each  epoch;  therefore,  the  control  gain  is  solved  for  once  and 
held  constant  throughout  the  simulation.  The  formation  slewing  Matlab  simulation  is 
given  in  Appendix  D. 

8.3.2  Formation  Slewing  Simulation  Results 

The  formation  slewing  simulation  provides  a  plethora  of  useful  results.  The  hub 
tracking  error,  estimation  error,  and  control  are  found  for  focal  lengths  of  both  0.5  and  4 
km.  Each  drones’  tracking  error,  estimation  error,  and  control  are  determined  as  well. 

For  conciseness,  only  the  results  for  two  drones  will  be  shown  (satellite  numbers  2  and 
31).  Two  slewing  angles  are  investigated,  90  degrees  and  30  degrees.  Figure  15  provides 
a  clear  image  of  the  entire  SI  formation  slewing  90  degrees,  with  a  0.5  km  focal  length. 
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y  (km) 

Figure  15.  SI  formation  before  and  after  90  degree  slew  (0.5  km  focal  length) 
The  magenta  circles  represent  drones  at  the  beginning  of  the  simulation,  and  the  red 
circles  represent  drones  at  the  end  of  the  simulation.  The  hub  is  the  black  asterisk  at  the 
origin.  The  upper-right  plot  illustrates  the  Golomb31  rectangle  formation  projected  into 
the  x-z  plane.  The  lower-left  plot  clearly  shows  the  drones  slewing  90  degrees  about  the 
hub-centered  x  axis. 

Figures  16-18  show  the  tracking  errors  for  the  hub  and  2  drones.  Although  the 
plots  are  specific  to  a  90  degree  slew  with  a  0.5  km  focal  length,  they  are  qualitatively 
representative  of  all  the  different  slew  angle  and  focal  length  scenarios. 
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Figure  18.  Drone  31  tracking  error  (90  degree  slew  and  0.5  km  focal  length) 

For  both  90  degree  and  30  degree  slews  with  either  a  0.5  or  4  km  focal  length,  the  hub 
tracks  to  within  50  meters  of  its  reference  position  and  to  within  5  millimeters  per  second 
of  its  desired  velocity.  The  drones  all  track  to  within  3  meters  of  their  desired  reference 
positions  and  to  within  1  millimeter  per  second  of  their  desired  zero  velocities. 

Table  6  lists  the  total  position  tracking  error  after  one  day  for  the  different 


scenarios  when  noise  is  turned  off. 


Focal  length 
(km) 

Slew  angle 
(deg) 

Hub  position 
tracking  error 
(m) 

Drone  2  position 
tracking  error 
(m) 

Drone  31  position 
tracking  error  (m) 

0.5 

90 

8.33e-6 

3.075e-6 

3.075e-6 

30 

3.05e-6 

1.126e-6 

1.126e-6 

4 

90 

7.35e-5 

2.713e-5 

2.713e-5 

4 

30 

2.69e-5 

9.93e-6 

9.93e-6 

Tal 

t»le  6.  Formation  slewing  position  tracking  errors 

with  no  noise 

For  all  scenarios,  neglecting  noise,  the  velocity  tracking  error  is  essentially  zero  after  one 
day.  If  noise  is  turned  off,  the  tracking  errors  go  asymptotically  to  zero,  as  expected  with 
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a  linear  quadratic  regulator  control  strategy.  Clearly,  the  noise  and  estimation  errors  have 
a  significant  effect  on  the  tracking  errors. 

Figures  19-22  show  the  estimation  error  results  of  a  dozen  simulations  for  the  hub 
with  varying  slew  angles  and  focal  lengths. 


Figure  19.  Hub  estimation  error  for  0.5  km  focal  length  and  30  degree  slew 
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Figure  20.  Hub  estimation  error  for  0.5  km  focal  length  and  90  degree  slew 
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Figure  21.  Hub  estimation  error  for  4  km  focal  length  and  30  degree  slew 
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Figure  22.  Hub  estimation  error  for  4  km  focal  length  and  90  degree  slew 
The  red  lines  are  the  three-sigma  values  calculated  by  Equation  8.8-8. 13.  The  hub’s 
steady-state  position  three-sigma  values  are  about  50  meters,  and  the  steady-state  velocity 
three-sigma  values  are  about  1  millimeter  per  second  for  all  scenarios.  The  estimation 
errors  are  within  the  three-sigma  values  with  few  exceptions.  The  three-sigma  values 
change  for  each  simulation  (because  the  noise  is  random),  but  the  change  cannot  be  seen 
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for  the  hub  because  the  order  of  magnitude  of  change  is  much,  much  less  than  their 
overall  value. 


Figures  23-26  show  the  estimation  errors  for  a  dozen  simulations  for  the  first 
drone  satellite  (satellite  number  2). 


Figure  23.  Drone  2  estimation  errors  for  0.5  km  focal  length  and  30  degree  slew 
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Figure  24.  Drone  2  estimation  errors  for  0.5  km  focal  length  and  90  degree  slew 
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Figure  25.  Drone  2  estimation  errors  for  4  km  focal  length  and  30  degree  slew 
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Figure  26.  Drone  2  estimation  errors  for  4  km  focal  length  and  90  degree  slew 
Drone  2’s  estimation  errors  are,  with  few  exceptions,  within  the  three-sigma  values  for 
all  scenarios.  The  three-sigma  value  change  from  one  simulation  to  another  can  be  seen 
in  the  drone  estimation  error  plots.  The  range  from  the  hub  to  the  drone  is  either  0.5  or  4 
km,  whereas  the  range  from  the  hub  to  the  earth  is  about  1.5  million  km. 

Figures  27-30  show  the  estimation  errors  of  a  dozen  simulations  for  the  last  drone 


(satellite  number  31). 
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Figure  28.  Drone  31  estimation  errors  for  0.5  km  focal  length  and  90  degree  slew 
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Figure  29.  Drone  31  estimation  errors  for  4  km  focal  length  and  30  degree  slew 


Figure  30.  Drone  31  estimation  errors  for  4  km  focal  length  and  90  degree  slew 
Drone  31’s  estimation  errors  are  within  the  three-sigma  values  with  few  exceptions.  The 
steady-state  position  three-sigma  values  are  less  than  0. 1  meters  for  each  drone  in  every 
scenario.  The  steady-state  velocity  three-sigma  values  are  less  than  1  micrometer  per 
second  for  each  drone  in  every  scenario. 

From  the  control  efforts,  the  directional  AV  can  be  determined  for  each  satellite: 
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1440  JT 

AV*  =  Yk| - 

*  S  1440 

(8.42) 

1440  .  Jrr 

av =YU  dT 

y  S1  1440 

(8.43) 

1440  J T 

C- 

(8.44) 

The  total  AV  is  found  from  the  directional  AV ’s  by  Equation  8.17.  The  differences 
between  Equations  8.14-8.16  and  Equations  8.42-8.44  are  the  number  of  maneuvers  over 
the  course  of  the  simulation  and  the  maneuver  interval.  The  formation  slewing 
simulation  runs  for  one  day,  with  one  maneuver  per  minute  (1440  maneuvers),  whereas 
the  Lissajous  orbit  simulation  runs  for  359  days  with  one  maneuver  per  day.  Table  7 
shows  the  average  AV ’s  for  a  dozen  simulations  for  the  various  scenarios. 

Hub  AV  (m/s)  Drone  2  AV  (m/s) 

1.0705  0.8271 

1.1355  0.9395 

1.2688  1.1189 

1.8570  2.1907 

Table  7.  Average  formation  slewing  AV ’s 
The  larger  the  angle  the  formation  slews  through,  the  more  AV  is  needed.  Also,  the 
larger  the  focal  length,  the  more  AV  required.  Table  8  shows  the  corresponding 
propellant  masses  needed  to  achieve  the  AV ’s  given  in  Table  7,  with  the  assumptions 
that  the  Isp  of  the  low-thrust  thrusters  is  10000  seconds,  the  initial  mass  of  the  hub  is  550 
kg,  and  the  initial  mass  of  each  drone  is  100  kg. 


Focal  Length  Slew  Angle 
(km)  (deg) 


0.5  30 


Drone  31  AV  (m/s) 


0.8307 


0.9587 


1.1315 


2.1932 
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Table  8.  Average  formation  slewing  propellant  masses 
When  the  noise  is  turned  off,  the  required  AV  and  propellant  mass  are  reduced 
significantly.  Table  9  shows  the  AV ’s  and  Table  10  shows  the  corresponding  propellant 
masses  when  noise  is  removed. 


Focal  Length 
(km) 

Slew  Angle 
(deg) 

Hub  AV  (m/s) 

Drone  2  AV  (m/s) 

Drone  31  AV  (m/s) 

0.5 

30 

0.0504 

0.0853 

0.0998 

0.5 

90 

0.1581 

0.2150 

0.2315 

4 

30 

0.4420 

0.5896 

0.6441 

4 

90 

1.3945 

1.9446 

1.9469 

Table  9.  Required  AV ’s  for  formation  slewing  cases  without  noise 


Focal  Length 
(km) 

Slew  Angle 
(deg) 

Hub  mprop  (g) 

Drone  2  mprop  (g) 

Drone  31  mprop  (g) 

0.5 

30 

0.2826 

0.0870 

0.1017 

0.5 

90 

0.8864 

0.2192 

0.2360 

4 

30 

2.4781 

0.6010 

0.6566 

4 

90 

7.8182 

1.9822 

1.9846 

Table  10.  Required  propellant  masses  i 

'or  formation  slewing 

cases  without  noise 
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8.4  Formation  Reorientation 


For  some  stars,  one  snapshot  from  the  SI  formation  will  not  provide  enough 
sampling  data  for  sufficient  resolution.  In  these  cases,  the  drones  must  rotate  90  degrees 
and  take  another  snapshot  In  Section  8.3. 1,  an  assumption  is  made  that  the  drones  are 
located  behind  the  hub,  in  the  negative  y  direction.  Following  this  assumption,  the  90 
degree  reorientation  will  be  a  rotation  about  the  y  axis.  The  aim  from  the  drones,  through 
the  hub,  to  the  desired  star,  is  maintained  with  such  a  reorientation. 

8.4.1  Formation  Reorientation  Simulation  Development 

The  simulation  development  for  the  formation  reorientation  is  nearly  identical  to 
the  formation  slewing  simulation  in  Section  8.3.1.  The  simulation  is  set  to  run  for  one 
day  with  a  maneuver  interval  of  one  minute.  The  place  in  the  Lissajous  orbit  is  variable. 
Generator  gives  the  dynamics  and  hub  reference  for  the  corresponding  place  in  the 
Lissajous  orbit.  The  continuous  state-weighting  matrix,  control-weighting  matrix,  and 
process  noise  strength  are  given  in  Section  8.3.1,  as  well  as  the  method  for  discretizing 
and  solving  for  the  control  gain.  The  control  gain  for  the  formation  reorientation 
simulation  is  identical  to  the  formation  slewing  simulation.  The  initial  conditions  are 
built  through  Equations  8.20-8.31. 

The  difference  between  the  formation  reorientation  and  formation  slewing 
simulations  is  the  desired  reference  states  of  the  drones.  The  formation  slewing  reference 
states  are  based  on  a  variable  slew  angle,  a ,  rotation  about  the  x  axis,  whereas  the 
formation  reorientation  reference  states  are  based  on  a  90  degree  rotation  about  the  y  axis. 
The  imaginary  central  drone  has  the  same  reference  conditions  as  initial  conditions. 
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Setting 


ft*  =-B* +ftcmt  (8.47) 

and 

(8.48) 

for  i  from  2  to  31,  gives  the  drones’  reference  positions.  The  rotation  is  90  degrees 
counterclockwise  about  the  y  axis  when  looking  back  from  the  hub.  Once  again,  zero 
relative  velocity  between  the  drones  and  the  hub  is  desired,  so  Equation  8.39  holds. 

The  control  and  estimation  algorithm  from  one  epoch  to  another  is  identical  to  the 
formation  slewing  algorithm.  The  formation  reorientation  Matlab  simulation  is  given  in 
Appendix  E. 

8.4.2  Formation  Reorientation  Simulation  Results 

The  formation  reorientation  simulation  provides  similar  results  as  the  Lissajous 
orbit  and  formation  slewing  simulations.  The  tracking  errors,  estimation  errors,  and 
control  efforts  of  the  hub  and  drones  are  determined.  From  the  control  effort,  the 
required  AV  and  corresponding  propellant  mass  can  be  found.  Figure  31  shows  the  first 
four  drones  of  the  formation  before  and  after  reorientation  (with  a  0.5  km  focal  length). 
Only  four  drones  are  pictured  for  clarity. 
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Figure  31.  Formation  reorientation  (0.5  km  focal  length) 

The  scale  is  in  kilometers  in  all  four  pictures.  The  plot  in  the  upper-right  shows  the 
projection  of  the  formation  in  the  x-z  plane.  The  magenta  circles  are  the  drones  at  their 
initial  conditions,  and  the  red  circles  are  the  drones  after  the  simulation.  The  formation 
appears  to  have  rotated  clockwise  90  degrees  about  the  y  axis,  but  the  hub  is  into  the 
page,  so  the  rotation  is  counterclockwise  when  viewed  from  the  hub. 

Figures  32-34  show  the  tracking  errors  for  the  hub,  the  first  drone  (#2),  and  the 
last  drone  (#31)  for  the  0.5  km  focal  length  case. 
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Figure  32.  Hub  tracking  error  (0.5  km  focal  length) 

The  hub  tracks  to  within  40  meters  of  its  reference  position,  and  to  within  8  millimeters 
per  second  of  its  reference  velocity  for  both  focal  lengths. 
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Figure  33.  Drone  2  tracking  error  (0.5  km  focal  length) 
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Figure  34.  Drone  31  tracking  error  (0.5  km  focal  length) 

Every  drone  tracks  to  within  4  meters  of  its  reference  position,  and  to  within  1 .5 
millimeters  per  second  of  its  desired  zero  velocity  for  both  focal  lengths. 

When  noise  is  turned  off,  the  satellites  track  much  better  than  when  the  noise  is 
included.  Table  1 1  shows  the  total  position  tracking  errors  for  the  various  satellites  after 


one  day,  when  noise  is  eliminated. 


Focal  length  (km) 

Hub  position 
tracking  error  (m) 

Drone  2  position 
tracking  error  (m) 

Drone  31  position 
tracking  error  (m) 

0.5 

2. 147e-6 

0.792e-6 

0.793e-6 

4 

2.146e-6 

0.792e-6 

0.793e-6 

Table  1 1 .  Formation  reorientation  position  tracking  errors  with  no  noise 


The  tracking  errors  go  asymptotically  to  zero,  and  the  velocity  tracking  errors  are 
essentially  zero  at  the  end  of  a  day.  Just  as  with  the  formation  slewing  simulation,  the 
noise,  and  in  turn  the  estimation  errors,  are  the  largest  reason  for  imperfect  tracking. 

Figures  35  and  36  show  the  hub  estimation  errors  for  a  dozen  simulations  with  the 
different  focal  lengths. 
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Figure  36.  Hub  estimation  errors  (4  km  focal  length) 

The  steady-state  x  three-sigma  value  is  about  than  30  meters.  The  steady-state  y  and  z 
three-sigma  value  is  about  50  meters.  The  steady-state  velocity  three-sigma  values  are 
about  1  millimeter  per  second.  The  hub  estimation  errors  stay  within  the  three-sigma 


values  except  for  rare  occasions. 


Figures  37  and  38  show  the  estimation  errors  of  a  dozen  simulations  for  the  first 
drone  (satellite  #2),  and  Figures  39  and  40  show  the  estimation  errors  of  a  dozen 
simulations  for  the  last  drone  (satellite  #31). 
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Figure  39.  Drone  31  estimation  errors  (0.5  km  focal  length) 
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Figure  40.  Drone  31  estimation  errors  (4  km  focal  length) 

For  any  drone  and  either  focal  length,  the  steady-state  position  three-sigma  values  are  less 
than  0.1  meters,  and  the  steady-state  velocity  three-sigma  values  are  less  than  le-6  meters 
per  second.  Also,  the  estimation  errors  stay  within  the  three-sigma  values  with  rare 


exceptions. 
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From  Equations  8.42-8.44  and  Equation  8. 17,  the  AV  can  be  determined,  for 
each  satellite,  from  the  control  effort  required  to  reorient  the  formation.  Table  12  gives 


the  average  AV ’s  from  a  dozen  simulations  to  reorient  the  formation. 


Focal  Length  (km) 

Hub  AV  (m/s) 

Drone  2  AV  (m/s) 

Drone  31  AV  (m/s) 

0.5 

1.0126 

0.8421 

0.8095 

4 

1.0133 

0.8496 

0.8190 

1 

fable  12.  Average  formation  reorientation  AV ’s 

The  focal  length  has  no  discernible  effect  on  the  AV  needed  to  reorient  the  formation. 
This  makes  sense  because  the  rotation  is  about  the  y  axis,  and  the  focal  length  is  assumed 
to  be  the  measurement  along  the  y  axis  from  the  hub  to  the  drones.  Table  13  gives  the 
propellant  masses  that  correspond  to  Table  12.  Once  again,  the  hub  is  assumed  to  have 
an  initial  mass  of  550  kg,  the  drones  have  initial  masses  of  100  kg,  and  the  Isp’s  of  the 
ideal  thrusters  are  10000  seconds. 


Focal  Length  (km) 

Hub  mprop  (g) 

Drone  2  mprop  (g) 

Drone  31  mprop  (g) 

0.5 

5.6771 

0.8584 

0.8252 

4 

5.6811 

0.8661 

0.8349 

Table  ' 

3.  Average  formation  reorientation  propellant  masses 

Without  noise,  the  AV  needed  to  reorient  the  formation  is  much  less,  as  shown  in 

Table  14. 

Focal  Length  (km) 

Hub  AV  (m/s) 

Drone  2  AV  (m/s) 

Drone  31  AV  (m/s) 

0.5 

0.0408 

0.1529 

0.1496 

4 

0.0408 

0.1529 

0.1495 

Table  14.  Required  AV  for  formation  reorientation  without  noise 
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Table  15  gives  the  propellant  masses  that  correspond  to  Table  14. 


Focal  Length  (km) 

Hub  mprop  (g) 

Drone  2  mprop  (g) 

Drone  31  mprop  (g) 

0.5 

0.2287 

0.1623 

0.1525 

4 

0.2287 

0.1623 

0.1524 

Table  15.  Required  propellant  masses  for  formation  reorientation  without  noise 
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9.  CONCLUSION 


In  Section  2,  the  equations  of  motion  for  the  circular  restricted-three  body 
problem  were  investigated,  and  the  libration  point  concept  was  addressed.  By  reducing 
the  circular  restricted  three-body  problem  to  collinear  libration  points,  and  linearizing  the 
equations  of  motion,  the  problem  was  simplified  enough  to  apply  it  to  a  linear-quadratic- 
regulator  control  scheme. 

In  Section  3,  the  optimal  control  was  found  with  the  linear-quadratic-regulator 
method.  By  making  a  judicious  infinite  horizon  assumption,  the  control  gain  becomes 
constant.  In  Section  4,  the  continuous  system  was  sampled  to  model  a  more  realistic 
discrete  system,  and  the  discrete  optimal  control  law  was  found.  In  Section  5,  the  system 
was  augmented  for  multiple  satellites,  and  the  fundamentals  of  formation  flying 
explained. 

In  Section  6,  an  observer  was  introduced  to  estimate  the  states  of  a  system  by 
using  incoming  measurements.  The  optimal  method  of  updating  the  estimates  was  given, 
to  reduce  the  unwanted  process  and  measurement  noise.  Hence,  the  Kalman  filter  was 
derived.  The  Kalman  filter  was  adapted  to  accommodate  nonlinear  measurements,  with 
special  emphasis  on  the  range,  azimuth,  and  elevation  measurements  common  to 
spaceflight. 

Section  7  gave  a  background  to  the  Stellar  Imager  mission,  and  its  performance 
and  design  requirements.  The  preliminary  formation,  orbit,  and  control  designs  were 
explored  in  detail.  Finally,  in  Section  8,  a  reference  orbit  and  associated  dynamics 
specific  to  the  sun-earth  L2  point  were  determined.  A  more  realistic  orbit  and  associated 
dynamics  were  found  with  the  use  of  Generator.  Three  simulations  were  developed  to 
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analyze  various  aspects  of  the  Stellar  Imager  mission.  Using  the  control  and  estimation 
methods  developed  throughout  the  thesis,  the  simulations  revealed  preliminary  tracking 
errors,  estimation  errors,  and  required  control  efforts  for  all  satellites  in  the  formation. 
From  the  control  efforts,  the  necessary  AV ’s  and  propellant  masses  for  orbit 
maintenance,  formation  slewing,  and  formation  reorientation  were  determined. 

The  control  strategy  and  Kalman  filter  provided  satisfactory  results.  The  hub 
satellite  tracks  to  its  reference  orbit  sufficiently  for  the  SI  mission  requirements.  The 
drone  satellites,  on  the  other  hand,  track  to  only  within  a  few  meters.  Without  noise, 
though,  the  drones  track  to  within  several  micrometers.  The  first  tier  of  the  proposed 
control  scheme  for  SI  requires  the  drones  to  track  within  centimeters.  This  could  be 
accomplished  with  better  sensors  to  lessen  the  effect  of  the  process  and  measurement 
noise.  Tuning  the  controller  and  varying  the  maneuver  intervals  should  provide 
additional  savings  as  well.  The  Kalman  filter  performed  such  that  the  estimation  errors 
were  for  the  most  part  within  the  three-sigma  values.  The  propellant  mass  and  AV 
results  provide  a  minimum  design  boundary  for  the  SI  mission.  Additional  propellant 
will  be  needed  to  perform  all  attitude  maneuvers,  tighter  control  requirement  adjustments, 
and  other  mission  functions. 

In  the  future,  many  research  opportunities  exist  for  formation  flying  satellite 
control,  libration  point  dynamics,  and  the  Stellar  Imager  and  other  similar  missions. 
Specifically  for  SI,  the  attitude  dynamics  and  control  problem  must  be  examined  and 
integrated  into  the  translation  control  problem  explored  in  this  thesis.  Also,  non-ideal 
thrusters  must  be  modeled  as  well.  In  simulation,  an  implicit  assumption  is  made  that  the 
thrusters  give  perfect,  impulsive  accelerations  determined  by  the  control  law. 
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Realistically,  thrusters  are  misaligned,  non-impulsive,  and  have  upper  and  lower  bounds 
on  force  outputs.  Collision  avoidance  is  another  problem.  In  the  simulations,  the  drones 
are  imaginary  points  that  may  pass  through  each  other.  In  actuality,  collisions  would 
cause  catastrophic  satellite  failures.  System  reliability  and  determination  must  be 
accomplished  as  well  How  a  drone  failure  is  detected  and  how  the  formation  responds  to 
such  failures  are  questions  that  must  be  addresed.  Nonlinear  control  and  estimation 
approaches  should  be  considered  also.  For  the  nanometer  level  preciseness  of  the  second 
and  third  control  tiers,  unique  control  strategies  and  algorithms  must  be  developed  and 
tested. 
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Appendix  A:  Generator-to-Matlab  Reference  Orbit  Conversion  Script 

getref.m 


disp ( ' loading  reference ' ) ; 

posbig=textread ( ' pos . txt '  ,  *  %s ' , ' delimiter ' , ' \n 1 ) ; 
%pos . txt  is  Generator  output 

velbig=textread ( ' vel . txt ' ,  '  %s ' , '  delimiter ' , ' \n  * ) ; 
%vel.txt  is  Generator  output 
for  i=l:359, 
for  j=l:3, 

p(  j ) =str2num(posbig{ (i-1) *3+j } ) ; 
v( j ) =str2num(velbig{ (i-1) *3+j } ) ; 

end 

v=v*  86.4; 

stater{i}=[p' ;v' ] ; 

end 

disp ( ' reference  loaded ' ) ; 
clear  posbig  velbig  p  v  i  j 
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Appendix  B:  Generator-to-Matlab  Dynamics  Conversion  Script 

getA.m 


disp ( *  loading  dynamics ‘ ) ; 

Abig=textread ( 1 Amatrices2 . txt  * , ' %s ' , ' delimiter ' , ' \n ' ) ; 
%Amatrices2 . txt  is  Generator  output 
for  i=l:359, 
for  j=l:36, 

a( j ) =str2num(Abig{ (i-1) *87+51+j } ) ; 

end 

Agen{i}=reshape(a# 6, 6) 
end 

disp ( ' dynamics  loaded ' ) ; 
clear  Abig  a  i  j 
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Appendix  C:  Lissajous  Orbit  Matlab  Simulation 
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%SYSTEM  DYNAMICS  AND  COST  WEIGHTING 
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pos { el } =s tatebarl { i } (6*el-5 : 6*el-3 ) ; 

range (el) = (transpose (pos {el} ) *pos{el} ) A .  5; 

elev(el) =asin(-pos{el} (3) /range(el) ) ; 
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plot (statem( end, 6*i-5) , statem(end, 6*i-3 ) , ' or ’ 
hold  on; 

plot (statem(l, 6*i-5) , statem(l, 6*i-3 )  , '  om' ) ; 
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subplot (3 ,2 ,4) , plot ( t , (statedifm ( : , 5) -stateestm( : ,5) ) *1000/86400, t, 3 *sigmaydotm( :,1) *10 00/8 64 00, 
3*sigmaydotm{ : , 1) *1000/86400, 'r ' ) ; 
ylabel('ydot  estimation  error  (m/s) ’ ) ;grid  on; 

subplot (3 ,2,6) , plot ( t , ( statedifm ( : , 6) -stateestm( : , 6) ) *1000/86400 , t, 3 *sigmazdotm( :,1) *1000/86400, 
3*sigmazdotm( :,l)*1000/86400, *  r  ’  )  ; 


subplot (3,2,6) ,plot (t, (statedifm( : , 6*i) - 

stateestm( : ,  6*i) ) *1000/86400, t, 3 *sigmazdotm( :,i) *1000/8 64 00, 'r'  ,  t, -3*sigmazdotm( :,i) *1000/864 00 
ylabel(’zdot  estimation  error  (m/s) ') ;xlabel (' time ') ;grid  on; 


Title:  Formation  Flying  Satellite  Control  Around  the  L2  Sun-Earth  Libration  Point 
Author:  Nicholas  H.  Hamilton 

Rank/Service:  Second  Lieutenant,  United  States  Air  Force 

Year:  2001 

Number  of  pages:  165 

Degree:  Master  of  Science 

Institution:  The  School  of  Engineering  and  Applied  Science  of  the  George  Washington 
University 


Abstract: 

A  growing  interest  in  formation  flying  satellites  demands  development  and 
analysis  of  control  and  estimation  algorithms  for  station-keeping  and  formation 
maneuvering.  This  thesis  discusses  the  development  of  a  discrete  linear-quadratic- 
regulator  control  algorithm  for  formations  in  the  vicinity  of  the  L2  sun-earth  libration 
point.  The  development  of  an  appropriate  Kalman  filter  is  included  as  well.  Simulations 
are  created  for  the  analysis  of  the  station-keeping  and  various  formation  maneuvers  of  the 
Stellar  Imager  mission.  The  simulations  provide  tracking  error,  estimation  error,  and 
control  effort  results.  From  the  control  effort,  useful  design  parameters  such  as  AV  and 
propellant  mass  are  determined.  For  formation  maneuvering,  the  drone  spacecraft  track 
to  within  4  meters  of  their  desired  position  and  within  1.5  millimeters  per  second  of  then- 
desired  zero  velocity.  The  filter,  with  few  exceptions,  keeps  the  estimation  errors  within 
their  three-sigma  values.  Without  noise,  the  controller  performs  extremely  well,  with  the 
drones  tracking  to  within  several  micrometers.  Each  drone  uses  around  1  to  2  grams  of 
propellant  per  maneuver,  depending  on  the  circumstances. 
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